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§ 1. Introduction. — § 1*0. The object of this paper is to develop methods whereby 
the differential equations of physics may be applied more freely than hitherto in the 
approximate form of difference equations to problems concerning irregular bodies. 

Though very different in method, it is in purpose a continuation of a former paper 
by the author, on a '^ Freehand Graphic Way of Determining Stream Lines and 
Equipotentials " (^ Phil. Mag.,' February, 1908; also ' Proc. Physical Soc.,' London, 
vol. xxi.). And all that was there said, as to the need for new methods, may be taken 
to apply here also. In brief, analytical methods are the foundation of the whole 
subject, and in practice they are the most accurate when they will work, but in the 
integration of partial equations, with reference to irregular-shaped boundaries, their 
field of application is very limited. 

Both for engineering and for many of the less exact sciences, such as biology, there 
is a demand for rapid methods, easy to be understood and applicable to unusual 
equations and irregular bodies. If they can be accurate, so much the better ; but 
1 per cent, would suffice for many purposes. It is hoped that the methods put 
forward in this paper will help to supply this demand. 

The equations considered in any detail are only a few of the commoner ones 
occurring in physical mathematics, namely : — Laplace's equation V^^ = ; the 
oscillation equations (V^+F)^ = and (V^—F)^ = 0; and the equation V*^ = 0. 
But the methods employed are not limited to these equations. 

The Number of Independent Variables. — In the examples treated in the paper this 
never exceeds two. The extension to three variables is, however, perfectly obvious. 
One has only to let the third variable be represented by the number of the page of a 
book of tracing paper. The operators are extended quite simply, and the same 
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methods of successive approximation apply. But, of course, the labour would be 
greatly increased. 

§ 1*1. Finite differences have, in themselves, but little importance to the student of 
matter and ether. They are here regarded simply as a makeshift for infinitesimals ; 
and the understanding is always that we will eventually make the differences so 
small that the errors due to their finite size will be less than the errors of experiment 
or practical working, and may therefore be disregarded. That it is possible to make 
them small enough without much labour is illustrated by examples given hereafter. 

In consequence of this point of view, the notation employed for finite differences is 
very similar to that for infinitesimal differences. Thus d and 3 are differential 
operators, while 8 and 'S are the corresponding finite difference operators. The oft- 

gS g2 02 ^ 

occurring symbol V^ = ^-^ + ^-^ + —^ is represented in finite differences by ^ ^. 

The differences employed are '' central differences," that is to say, they are 
considered as existing at the centre of the group of co-ordinate points from which 
they are derived. In this respect the notation differs from that used and defended 
by Boole (^ Calculus of Finite Differences,' Art, 14), in which A\. is considered to 
exist at one end of the set of n+1 quantities which contribute to its value. The 
differencing operator S and Sheppaed's"^ averager /x are defined by 

S/H =/(^+P)"~"/(^™P) • • (1). 

m/(^) =/(^ + P) + /(^-P) (2), 

where li is the co-ordinate difference. Sheppaed shows that /x, S, and d combine with 
one another according to the ordinary rules of algebra. 

In this paper the differential coefficient d''f{x)\dx''' will be approximately represented 
when n is even by h~'' .h""/ {x) at the tabular points, and by /^~^/xSy(a^) half-way 
between the tabular points. That is to say, these difference ratios are taken in place 
of the differential coefficients, and the error caused by so doing is left for consideration 
until after the difference equation has been solved. When n is odd, the symbolic 
expressions given above for d''f(x)/dx'' at tabular and at half-way points are simply 
interchanged. The representation is closer when the averager jx need not be intro- 
duced. Partial differential coefficients are represented by the difference ratios found by 
performing the above operations with [jl and 8 for each independent variable in turn. 
It will be convenient to have the representation of some of the commonest differential 
coefficients set forth explicitly. 

Let ^ be a function of x and t/, and let lines be ruled on the plane xy parallel to 
the axes at equal distances /^, of x and y, so as to divide the surface into a number of 
equal squares each of side h units. Let the arithmetical value of cj) at the centre 
point of each square be written down in the square, forming a table of double entry. 

* W. F. Sheppard, " Central-Difference Formulae," ^ Proc. Lond. Math, Soc.,' vol. xxxi., p. 460. 
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Let (01), (23)5 &C.5 for examplej represent these numerical values at the co-ordinate 
points X = 0, y = --h; x = 2h, y = --3/^5 &c., reckoned from the centre of an 
arbitrarily chosen square as origin. Then at any tabular point 0, 



dx 
dx^ 

3^ 

dx^ 

3^ 

dx^ 

dxdy 
dx^ dy^ 



will be represented by 



J? 



?? 



?? 



?5 



?? 



?? 



J? 



?? 



?? 



J? 



^ = ^{(lO)-2(00) + (IO)} 
^.= ^.{(20)-4(l0) + 6(00)-4(T0) + (20)} 
.. §+^= ^/<^ = i{{10) + (0l) + (l0) + (0l)-4(00)}. 



Thus h^^^f is the sum of the four nearest neighbours minus four times the value 
at the point considered. 

^1'^ = ( DZi + 2 5::j^:::2+ ^4)<^ ^i^l be represented by (^-^ + 2 ^^-;^^^ + ~)(i = ^i^(^ 



^dx^ dx^ dy^ dy^ 



\^x' ^x'^y' % 



and ^iV = |j[20(00)-8{(10) + (01) + (T0) + {01)} + 2{{ll) + (ll) + {ll) + (ll) 



Halfway between two tabular points, say at |-, 



+ {(20) + (02) + (20) + (02)}]. 



^ will be represented by ^^ = y{(lO)~(00)}. 



dx 



"^x h 



In. the centre of four tabular values, e.g.^ at ^^ ^ 

^ will be represented by ^ = l{(li) + (00)-(0l)-(10)}. 

A point at which the difference equation obtaining throughout the body has to be 
satisfied will be called a body-point. There must be enough known values of the 
integral <j) on the boundary side of any body-point to make the said difference 
equation completely determinate. Thus for ^ ^ there must be at least one layer of 
points with known values of <^ on all sides of any body-point, for ^ ^ at least two 
layers. It will be seen in § 4 that at re-entrant angles a point may have sufficient 
known values outside it and yet not be a body-point, because the body-equation is 
not to be satisfied there. A point at which the body-equation is not satisfied, but at 
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which there is a value of cj) which enters into the system of body-equations by way. of 
the boundary-conditions will be called a boundary -point 

The values of any function of position at these two classes of points will be 
distinguished as body values and boundary values^ or synonymously as body-numbers 
and boundary-numbers. 

Problems are divided into two maiii classes according as the integral can or cannot 
be stepped out from a part of the boundary. They are discussed in § 2 and § 3 
respectively. 

§ 1*2. Errors due to Finite Differences, — Having solved an equation using the 
simple expressions of § 1*1 for the differential coefficients, it remains to enquire how 
much in error the integral may be. A rule of apparently universal application is to 
take smaller co-ordinate differences and repeat the integration ; and, if necessary, 
extrapolate in the manner explained below. 

It is known"^ that when central differences are used, the expansions of the differential 
coefficients of a function in terms of its differences contain only alternate powers of 
the co-ordinate difference h. The same is true for partial differential coefficients and 
for products of differential coefficients. Consequently the error of the representation 
of any differential expression by central differences is of the form h^V^ (^j Vr ^) 
+ /^^F4 {x^ y^ z)-i- terms in higher powers of h^^ where Fg, F4, &c., are independent of JiA 

Next, as to the error of the finite-difference-integral ^. This is the infinitesimal 
integral of a differential equation having the error h^l^^ (^j y^ ^) + ^^^^4 (^? 2/> ^) + ? &c. 
Let F be the integral of the correct differential equation. Then, if we write 

<l> (^, y, z) = F (x, y, z)+m^i {x, y, z)+mP\lj2 (x, y, z) + terms 

in higher powers of m, it follows that a differential expression of any order and degree 
for (jy differs from the corresponding one for f by 

m X (a function of the differential coefficients of F and of t//i)+ terms in m^, m^, &c., 

provided only that m is independent of the co-ordinates. Now, identifying m with h^, 
it follows that : the errors of the integral and of any differential expressions derived 
from it, due to using the simple central differences q/' § 1*1 instead of differential 
coefficients^ are of the form 

h% {x, y, z)+h% {x, y, z) + h% {x, y,z) + , &c. 

Consequently, if the equation be integrated for several different values of h, 
extrapolation on the supposition that the error is of this form will give numbers very 
close to the infinitesimal integral. When h is small enough the error is simply 

^ W. F. Sheppard, '« Central-Difference Formulae," 'Proc. Lond. Math. Soc.,' vol. xxxi. (1899). 

[t Mte added January 21, 1910. — It is assumed that the co-ordinate axes in the tables which are 
compared are parallel, for the error at a ixed point and for a fixed value of h may depend on the direction 
of the axes.] 
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proportional to h^. Peculiarities present themselves on the boundary, but it is easy 
to see that errors will be of the form h^f2 + hy4: + , &c., provided that in passing from 
one table to another each part is either kept infinitesimally correct, or else is worked 
by differences whose size in the one table bears a constant ratio to that in the other. 

An extrapolation can only be made where the tabular points of the several tables 
coincide with one another. It is conceivable that in the future some method 
will be found of defining a continuous function in terms of the discrete body and 
boundary values, so that this continuous function shall have an error of the form 
hY2{^iyi^)-^hy^{x,y, z) + , &c., everywhere. Extrapolation would then be possible 
everywhere. 

An excellent illustration is afforded by Lord Rayleigh's account of the vibration of 
a stretched string of beads {' Sound,' vol. I., § 121). He gives the frequency of the 
fundamental for the same mass per unit length concentrated in various numbers 
of beads. This is reproduced below in the table. The co-ordinate difference h is 
inversely as one plus the number of beads, not counting beads at the fixed ends. 



Number of free beads + one . . 


2 


3 


4 


5 


10 


20 


40 


00 


Ratio of frequency to that of"! 
continuous string .... J 


•9003 


•9549 


•9745 


•9836 


•9959 


•9990 


•9997 


unity 


Error in representation of con-1 
tinuous string by string of > 
beads 


•0997 


•0451 


•0255 


•0164 


•0041 


•0010 
•4111 


•0003 


•0000 


Ratio of error to square ofl 
co-ordinate difference x a > 
constant . . . . . . .J 


•3988 


•4059 


•4080 


•4091 


•4107 


•4112 


•4112 



The degree of constancy of the last line shows that if we found the frequency for 
one bead and for three, then extrapolation, on the assumption that the error is 
proportional to h^, would give us the frequency for the continuous string to about one 
part in 1000 ; which is as near as we could get by twenty beads and no extrapolation. 
While extrapolation from the exact solutions for four beads and for nine would leave 
an error of only one in 50,000. Other examples of extrapolation will be found in § 3*1. 

§ 2. Procediire %vhen the Conditions alloiv the Integral to he Marehed out fro'fn a 
Part of the Boundary. 

§ 2 ' 0. Histmical. — Step-by-step arithmetical methods of solving ordinary difference equations have long 
been employed for the calculation of interest and annuities. Kecently their application to differential 
equations has been very greatly improved by the introduction of rules allied to those for approximate 
quadrature. The papers referred to are : — 

* * 

KuNGE, " Tiber die numerische Aufiosung von Differentialgleichungen," ' Math. Ann.,' Bd. 46. 
Leipzig, 1895. 



312 ME. L. F. EICHARDSON: APPEOXIMATE AEITHMETIOAL SOLUTION 

"W. F. Sheppard, ** a Method for Extending the Accuracy of Mathematical Tables," ' Proc. Lond. 

Math. Soc./ XXXI. 
Karl Heun, "Neue methode zur approximativen Integration der DifFerentialgleichungen einer 

unabhangigen Verander lichen," * Zeitschrif t Math. u. Phys.,' No. 45, 1900. 
WiLHELM KUTTA, " Beitrag zur naherungsweisen Integration totaler Differ en tialgleichen," ' Zeit- 

schrift Math. u. Phys.,^ No. 46, 1901. 

Further Eichard Ganz, in a paper " Uber die numerische Auflosung von partiellen DifFerential- 
gleichungen," ' Zeitschrift Math. u. Phys.,' No. 48, 1903, has extended the methods of Eunge, Heun, and 
KuTTA to partial equations of the type considered in this section. Those of the first order he turns into 



^ = 'F (X, 7/, 4>, ^ 

ox \ cy 



and starting from a boundary, where </> is a given function of z/, he expands this function as a power series 
of y and integrates step by step in the x direction. The results he gives are of remarkable accuracy. 
It is less accurate, but simpler, to work entirely by arithmetic in the manner illustrated in § 2*2 below. 

§ 2'1.'^ A simple Process and its Possibility, — For ordinary equations the necessary 
and sufficient condition is that, for an n^^ order equation, the integral and all its first 
n—l differential coefficients should be given at the boundary. This is almost obvious 
at first sight. The complications in the following arise entirely from having to attend 
to the correct centering of the differences—an important thing in practice. Let the 

equation be ^ ^ f / /. # ^ d''~'^ \ 

Then if all the quantities of which f is a function are given at a^ = Xq^ we can 
calculate d^'cjy/dx'' at ^o- Now, representing differentials by simple central differences, 
draw up a table in columns. The subscripts denote distance from Xq, 

xo + h 4>h {^^)h 

Xo + fh m)Bi2h ' 

Xo-\-2h (^2h 

Each difference is centred at values of x halfway between those for the difference 
of next lower order. Then, beginning with S^'i^, each difference is added to the one of 
next lower order and the sum written down in the column of the lower order one, a 
step A after it. This process gives a table with a diagonal boundary. The next 
value of S""^ is found from the difference equation and the process is repeated. The 
difference equation is satisfied at those values of x where the highest difference is 
tabulated. These processes involve moving certain differences through a step of ^h^ 
at the start and when satisfying the difference equation. Except for S"" and S"""^ this 
is done by (8^^)i/^^ = (8'*<^)^ + |-^ (8'*"^^^)o. Note that this is a central formula with 
respect to the step —^h to -{-^K For S"""^ and S"" the motion of^h is accomplished by 

^ Eevised April 20^ 1910 
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writing algebraic symbols for the unknown values of S'^c^ and S'*""^<^ and finding them 
from SS''~'^(f) — S""^, and from the given difference equation simultaneously. On the 
other hand if any of the first (n — l) diflPerential coefficients are unknown at x^ the 
summation from column to column across the table cannot be carried out, and so the 
step method is impossible. The conditions for stepwise integration of partial equations 
are not so simple and require investigation. 

§2*2. As an example of stepwise integration^ let us take the equation for diffusion 
in a tube, cr9^(^/3x^ = d(f)/dL In the first place, if <^ = F {x, t) is a solution of 
d^^/dx^ = d^/dt, then (j) = ¥ {ax, a^crt) is a solution of the given equation. So we 
need only concern ourselves with d^cfy/dx^ = d(f)fdt, and the results will apply to bodies 
of any linear dimensions and any uniform diffusivity. Let us suppose the boundary 
conditions are : ^ = when x = ±f for all values of i^ ; <j& = 1 for all values of cc when 
^ = 0; in fact the familiar case of a uniformly heated slab, the faces of which are 
suddenly cooled. The method of this example is so simple that it can hardly be 
novel. It is introduced to show how easy it sometimes is to obtain approximate 
integrals by arithmetic of equations usually treated by complex analysis. We draw 
up a table with a row for each 0*1 oi x and a column at every O'OOl of t, (The reason 
for making the time step small will appear later, § 3 •2*1.) The given boundary values 
are next inserted. 

Table I. 





t - 0. 




0-001. 


0-002. 


0-003. 


0-004. 


- 005. 


X = 0-5 


• 0000 




0-0000 


0-0000 


0-0000 


- 0000 


0-0000 


0-4 


1-0000 


a 


= 0-9090 


0-8356 


0-7714 


0-7209 


0-6729 


0-3 


1 • 0000 


h 


- 0-9959 


0-9834 


0-9695 


0-9492 


0-9329 


0-2 


1-0000 


c 


- 0-9998 


0-9993 


0-9968 


0-9945 


0-9887 


0-1 


1-0000 


d 


- 1-0000 


1 • 0000 


0-9999 


0-9994 


0-9990 


0-0 


1 - 0000 


e 


- I -0000 


1-0000 


1 • 0000 


1 - 0000 


0-9998 


-0-1 


1 • 0000 


d 


= 1-0000 


1 • 0000 


0-9999 


0-9994 


0-9990 



0-005 

correct by 

Fourier's 

method. 



Errors. 








1 




0000 
6828 
9545 
9980 
9996 
0000 
9996 



-0-0099 
-0-0216 
"0-0093 
-0-0006 
~ - 0002 
-0-0006 



In satisfying the equation we must be careful to equate values of 'S^<^/'Scc^ and 
^4>/^t, which are centered at the same point. This causes a little difficulty at starting. 
When t = 0*001 let the values of (f) be a, b, e, d, e, as indicated in Table I. Then if 
the difference equation be satisfied at ^ = 0*0005, it takes the form of 5 simultaneous 
equations involving a, b, c, d, e. Solving these equations, we find the numbers given 
in the column x = 0*001. Having got over this rather troublesome first step, we can 
find the rest much more simply by centering all differences on the columns t = 0*001, 
0*002, 0*003, &c., and deducing each number from the two preceding columns. The 
errors resulting from the above process may be found by comparison with the Fourier 



ni- 



solution ^= S ( — 1)^ ~e ""^'^'^ cos (m7r^). This series has been computed when 

m odd m; 

t ■=. 0*005, and the numbers so found are given in the table. It is seen that the step 
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method gives a fair general view of the temperature corresponding to a given time and 
position, and it can, of course, be employed with equal facility when the boundary 
temperatures vary in almost any assigned manner, including cases that would be 
difficult or impossible by the Fourier method. 

§ 3. It frequently happens that the integral ean only he determined with reference 
to the boundary as a tvhole, as, for example, in the calculation of the electrostatic 
potential at all points of a region when its value is given over the conducting 
boundaries. Here the differential equation is of the second order, and the first space 
rates of its integral are not given on the boundary, so that the step-method is 
inapplicable. The following § 3 contains an account of two methods for solving 
problems of the type indicated. 

§ 3'0. The Determinate Nature of the Pro&Zem.— Let there be n body-points and 
s boundary-points in the region considered. Then the differential equation, to be 
satisfied in the body, is approximately represented at any body-point by an algebraic 
equation connecting the body value there with the surrounding values. This algebraic 
equation will be of the first, second, or higher degree, according as the differential 
equation is of the first, second, or higher degree in the function of position and its 
difierentials. Forming this equation at every body-point, we have a system of 
n simultaneous integral equations between s + n unknowns. To make the problem 
determinate, the boundary conditions must therefore supply s independent relations, 
involving the boundary values. The rules governing the 
arrangement of these 5 boundary equations, so as best to Table 11. 

represent the given infinitesimal boundary conditions, have 
not yet been elaborated. In certain cases a choice of ways is 
open, as in the following example : Let the body equation be 
(ff^/dx^-\-d^/dy^) f ~ /), where p is a given function of position? ., 

and the boundary condition df/dn + ^f ^ 0, where ^ is an ^ 

arbitrary function of position on the boimdary. Let the body 

values of/ be denoted by \//i, \//2, ..., and the boundary values 

by /3i, ySa, &c. Then in the annexed Table II. we are at liberty to choose between 

two alternative approximations. For we may take as values of / on the boundary 

And as corresponding values of ^f/^n, 

A-V'lJ ^2-'^2, A~^2, A-^S. 

So that we have one boundary condition for each /3, Or else we may suppose the 
corner slightly bevelled, so that, while the above relations still hold for ^i and yS^, 
we now have at the corner f = i {^2 + ^(^2 + /Ss)}, ^f/^n = ^{/32 + /3s)-xfj2. Here 
^2+1^3 enters as one variable. Also /Sg and ySg only enter the body equations 
when combined in the form ySs + ySg. We no longer seek to determine ^82 and ^83 
separately, but merely their average, which is the value of/ just outside the corner. 
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so that there is still one boundary equation for each boundary variable which we seek 
to determine separately. We have just considered a blunted angle of 90 degrees. By 
examining the remaining five angles possible with two co-ordinates, namely, 45, 135, 
225, 270, 315 degrees, one may convince oneself that with a closed boundary it is 
-always possible to arrange to have just as many boundary equations as boundary 
unknowns. It is in some cases necessary to suppose the corner slightly bevelled 
and to replace certain (^)'s by their averages. The representation of the boundary 
condition when both f and dfjdn are given at each point will be considered in the 
theory of the dam. 

When the infinitesimal boundary conditions are such as to make the problem 
determinate, it will be assumed that we can, and therefore do, represent them by a 
set of boundary equations equal in number to the boundary unknowns. If any case 
be discovered in which this is impossible, it will be an exception to the rest of § 3. 

§ 3*1. The finite difference problem being thus made determinate the most direct 
way of finding the integral is to solve the n-\-s simultaneous algebraic equations for 
the n body and s boundary values of the integral. To take an example : — At 
one pair of opposite 

edges of a square </.= !, Table III. 

at the other pair ^ = 0. 
Inside 



0-5 



1 



1 



1 



edge 

of 

square. 



d'<i>idx'~\-dmdf = 







0-5 



a 



h 



c 



h 











d 



e 



d 







0-5 



/ 



0-5 



everywhere. Find ^ in- 
side the square. 

Now by symmetry 
the values of ^ on the 
diagonals will be every- 
where 0*5. In fact, we 
need only consider |- of 
the area 'of the square ; 
all the rest follows 
from it. Takinsr finite 
difierences, Table III. 
is drawn uj) with the 
given boundary and 
diagonal numbers in 
their proper places and 
a, b, e, d, e, f for the 

unknown body values. Now as the finite difierence expression for V^(j) given in § 1*1 
has to vanish at all body-points, we have a relation between each of the letters 
a to f and its four nearest neighboui^s. 

2 s' 2 







^ 



centre of 

0-5 
square. 



edge 
of square. 
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The solution of these six simultaneous equations was accomplished in an hour and 
gave the following results : — 



By finite differences 



By infinitesimals 



Errors due to finite differences 



• • 



a. 


k 


0-794 
0-798 
0-004 


d. 


e. 


0-533 


0-693 


0-112 


0-601 


0-6324 


0-700 


0-111 


0-604 


0-6354 


0-534 


0-007 


0-005 


0-003 


0-003 


0-001 



The numbers for infinitesimal differences were obtained from 



A Vl r 1 TYhTT 

d) =: ~ X/^i)2___ sech — cos mx cosh mz^ 

TT modd m 2 

the separate terms of which satisfy V^cf) = at all points, and ^ = when x = ±^7t, 
and by their addition make c^ = 1 when z = ±|-7r. Adding up the series at these six 
points took 3 hours. It is seen that the greatest error is 1 '4 per cent, of the range of 
potential between the side and diagonal. 

Further, if we take co-ordinate differences of twice this size, leaving only one 
unknown in the same position in the square as e here occupies, its value is easily found 
to be 0*6250. Extrapolating as in § 1*2 we find for infinitesimal differences at this 
point (f) = 0*6324 + i(0-6324~-0-6250) = 0*6348, and this is only ^^oth per cent, in 
error. To correct the other values we should have to halve the co-ordinate difference 
instead of doubling it, and this would require much more work. 

As a second example of the use of simultaneous integral equations, let us take the 
determination of the gravest period of vibration of a thin square plate with edges 
clamped in a plane. It is known (Love's ^ Elasticity,' ed. 1906, p. 469) that the 
displacement normal to the plate is of the form W cos (p^-fe), and W satisfies the 

equation f -^— +2 ^ 4- ^-^ ) W E V/W = c^W, where c'^ is given as a function ofp, 

the elastic constants, the thickness and the density. Now let us form a table such 
as (IV.) to represent W. In this table W is measured from the plane of the clamped 
boundary. The differential equation, when turned into finite differences, becomes a set 
of simultaneous equations, connecting each in turn of the unknown body values q, r, s, t 
with its twelve nearest neighbours. As the boundary numbers are all zero there are 
no constant terms in these equations, and they are only consistent when the 
determinant of the coefficients of q, r, s, t vanishes. There are a number of values 
of d^ which cause the determinant to vanish, and of these the smallest is that 
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belonging to the gravest mode of vibration. Let Wi ~ f {x, y) be the appropriate 
integral for a square of unit side. Then Wl =/(^/L, yjL) will bear a similar relation 

L^ W, L^' 



to a square of side L and to a new constant G^, Then Cl'^ = 



Wl L^ Wl 

so that L^Cl^ is independent of the length of side of the square. It will be interesting 
to notice how this constant con- 































s 



verges towards a limit as the 
number of co-ordinate differences 
in the side of the square is 
increased. 

The configuration chosen was 
one in which the sides of the 
square are at 45 degrees to the 
rows and columns of the table. 
This gives a sharper boundary 
than the parallel arrangement. 
The symmetry of the gravest mode 
reduces the number of unknowns. 
As well as the arrangement in 
Table IV., two smaller ones were 
also considered, namely, those 
formed by cutting off in turn 
its first and second outer layers. 
C^ was calculated in each, case 
from the determinant, by approximation where necessary. 

Collected results : — 



Side of square = L. 



Table IV. 
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Ci} for the gravest mode. 


L^Cl*. 


20 • 0000 
6-3058 
' 1-88843 


405 - 0000^ 
985-28 
1133-53 



We may attempt a closer approximation by assuming that the error in C^^ is 
inversely as the square of the number of co-ordinate difierences in the side of the 
square, in accordance with §1*2. From the 1st and 2nd values extrapolation gives 
L^Cl^ = 1311-69. And from the 2nd and 3rd L'(V = 1287*96. Or, if we assume 
that the error is of the form eJi^-VeJi^ and extrapolate from L^C^^ for the three values 
of h, we find U'Gj^ = 1282*62. The way in which the values of L^C^^ converge 



■^ Ee vised figures, Marcli, 1910. 
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indicates that this last one can hardly be more than 1 in 1000 in error. The precise 
analytical theory of the vibrations of a square plate with clamped edges does not 
appear to be given anywhere, but we may obtain an upper limit to L*Cj,* by the 
method in Rayleigh's '' Sound," § 89. Assuming W = {a^—x^f {¥—y^f^^ we have 



r+b 



\y'WJdxdy \ 



+a 



C*> \ymjdxdy \ WMa;%= (31-5000) (a-*+6"') + {18-0000)(a&) 



—a 



when a = h this becomes 1296*0/(2^)"^, which is greater than 1283/(2a)^ obtained by 
finite differences. 

§ 3 '2. Successive Approximation to the Integrals, --^id^\Aiig illustrated the use of 
simultaneous integral equations, let us pass on to methods which have this property 
in common : that starting from a table of numbers, correct at the boundary, but 
otherwise merely as near as one can guess, one proceeds by definite methods to 
modify this table and thereby to cause it to approach without limit towards the true 
finite-difference integral. 

Conditions, — ^The following methods of approximation have up to the present 
been applied only to a limited class of equations satisfying the conditions given 
below. 

Let /be an arbitrary function of position having n body values i|fi, t/i^, ..., ^„, and s 
boundary values ^i, ^2^ •••j A- -Let the differential equation to be solved be J)<^ = 0, 
where J) is a differential operator, together with such boundary conditions as 
make the problem determinate. Let D be approximately represented by the finite 
difference operator 2)^ so that the body equations are 

3)'i/^i = 0, 3)%=0, ..., 3)V^^ = . . . . . '. (1). 

Then in order that the following justification of the approximation method may apply, 
it will be shown in the Appendix that 5)' and the boundary-equations must be linear ; 
and the body- and boundary-equations must be the condition that a certain positive 
homogeneous quadratic function V of t//i, -^2, -*-, ^n is a complete ininimum. Also 
y, though otherwise arbitrary, is limited on the boundary to be the difference of two 
functions of position both of which satisfy the said boundary- conditions. Under 

A: = n 

these circumstances / can be expressed in the form /— 't A^P^, where P^ is an 

A; = 1 

integral of ( + 3)'— X/) (^ = 0, which satisfies the same boundary-conditions as/, 
X/ being a positive constant, the sign before 2)^ being the same as that of xfif, in 
3)^t/;;fc. The proof of this fact and of various other properties of the (P)'s will be 
deferred to the Appendix. The (P)'s may be called the principal or normal modes 
of vibration of the system. V is analogous to potential energy. Some of the 

♦ I am indebted to Prof. A. E. H. Love for pointing out this method and for giving me the numerical 
result for a square. 
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commoner possible forms of 3)^ will be found in the table on p. 354. As well as the 
equation 2)^^ = 0, the equation S^^ = {a given function of the co-ordinates), and 
the equation (3)'— X^) <^ = may be treated by these successive-approximation- 
methods. They will be discussed in order. 

§ 3*2'1. The Equation 5)^^ = 0. — The approximation process proceeds as follows. 
Let j>y, be the correct finite-difference integral. Let <^i be a function (that is a table 
of numbers) satisfying the correct boundary-conditions, but arbitrary as to its body 
values. Next calculate the body- values of j>2 by means of 

^2 ^^ ^1 — %~ 2) <^i ......... . (1) 

where a^ is a number to be fixed ; and fill in such boundary-values of <^2 ^s will satisfy 
the same boundary-conditions as ^^. The succeeding steps are each of the form 

for the body values, and by choosing the boundary values (f)„^+i is made to satisfy the 
correct boundary condition. These are matters of simple arithmetic. It will be 
shown that by the judicious choice of %, a^, ..., a^ it is possible to make ^t+i nearer to 
(j)^ than ^1 was. For since S' is linear and S^<^m = we have from (2) 

Now it is shown in the Appendix that ^,n—^u may be expanded in a series of integrals of 

J^Uu ^j m^ =: 2yx\.-^.ir^ iOu 

Then by (4) J)^ (c^^ - c^„) = + S A,X/P,. 

And therefore by (3) ^2~^t. = ^AJ l--^j P^. (6). 

Proceeding in the same manner after t operations we arrive at 

A measure of the. deviation of two functions from one another which is used in the 
theory of Least Squares is the sum of the weighted squares of their differences. On 
the same principle let us measure the error of ^^+i by 

-I^^+i ^^ ^ (y*^+i — Yu/ ^ ^ ......... (8), 

where S stands for a summation over the body points and I is a certain one-signed 
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function of position. (See Appendix.) Squaring both sides of (7), multiplying by I 
and making the summation, we have by a property of the (P)'s (see Appendix, 
equations (31) and (32)) 



bi (<p^ + l <pu} — SA;;. 



«l/\ «2/ \ «« 



(9). 



Now it has been found that by a judicious choice of %, ag? -m ^t? the quantity 

may be made small for all possible values of X/. 



^1/ \ ^2/ \ ^t/J 



(Thus fig. 1 shows this done for a set of seven (a)'s. This graph was arrived at by 
trial.) The error E^+i of ^^+i may therefore be made small in comparison with that 
of (j)i. This is possible because the values of X/ lie in a finite range ; corresponding 
to the fact that there are only a finite number of terms in the series SA^^P/^. A similar 
process will not work with the infinite series of sines, Bessel functions and other 
infinitesimal integrals of (5)— X^) P = 0. In choosing a^, %, ..., a^, a diagram of the 
kind shown in figs. 1 and 2 is a great help. In this we take for the abscissa a 
variable X^ which takes in turn the values Xi^, Xg^, ..., X/, and as ordinate we consider 

ci> = (l— — j(l )x ... x( 1 — — ). ...... (10). 

The value of o> at X^ = X/ is the ratio of the amplitude of the vibration P/, in the final 
approximation ^^^i to its amplitude in the initial guess. The individual factors 
(l—k^/a^) in 0) represent straight lines, all cutting the vertical axis at 6) = 1 and the 
horizontal axis at the points X^ = a^, otg, . . . , a^. By bringing any adjacent pair a^ and 
a^+i closer together the values of a> corresponding to the range of X^ between a,, and 
a^+i are diminished, provided that the other (c6)'s remain fixed. It follows that by 
judiciously spacing the (flt)'s along the horizontal axis and by taking a sufiicient 
number of such points (that is of approximations) the successive maxima and minima 
of CO can be made all less, in absolute value, than any finite quantity e however small. 
When this is done S (^^+1— <^i*)^I being equal by (9) to SA/cy/ where co^ is the value of 
a> at X^ = X/ must be less than e^SA/. That is to say, the ratio of the error of the last 
approximation (f>t+i to that of the initial guess ^, being {SI(^^+i — ^,,)7SI(^l~-^^,)^}^^^ 
is less than €, and € can be made very small. 

A knowledge of X/, X^^ . . . , X/ is not necessary, but it is necessary to know the 
limits within which they range, or limits enclosing these. For the lower limit we 
require an estimate of Xi^, which in the dynamical application is the square of the 
frequency of the gravest mode of vibration multiplied by a constant depending on 
density and elasticity. It is usually sufficiently close to take some boundary such as 
a rectangle or a sector of a circle, for which the frequency is known, and which fits in 
a rough way the irregular boundary under consideration. 
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Nextj to find the upper limits it is always possible to proceed as follows : — Make a 
guess at A^P„, that is to say^ write down as sharply oscillatory a set of numbers as 
possible for the body values, such for example as Table V. 
when extended to fill the region, and add to it such boundary Iable V. 

values as will satisfy the boundary equations. Call this guess 
Xi. Then xi = A^P«+ A„_iP„_.i + , &c. Now operate on xi 
with ^^ many times in succession, and after each operation ^ 
readjust the boundary values so as to satisfy the boundary 
equations. By this process the coefficient of P^ is increased -^ -- 1 4-1 

relatively to the coefficients of the other (P)'s in the 

expansion of the resulting table, because X/ is the greatest of the X^. So that the 
table approaches a multiple of P.^. From P« it is easy to find X/. In any case a 
rough approximation to X/ suffices. 

When the boundary values A • • • A vanish this labour is unnecessary, for then X/ 
cannot exceed (see Appendix) the greatest value of X^ pertaining to an integral of 
( + S)^--X^) ^ = with the given size of co-ordinate differences, and with ^ vanishing at 
infinity ; and this value of X^ depends only on the form of S)^ and may be calculated 
once for all. It will be denoted by X^l Thus when 2)^ = ^^/^x^+^^/^y^ considerations 
of symmetry show that the most oscillatory integral is Table V. extended similarly 
in all directions, and from this we find X^^ = 4:/{^xy -h i/{^yy. 

Having thus found limits between which X/, Xa^, ..., X/ must lie, it remains to 
choose the (a)'s so as to make o) small for all value of X^ in this range. In practice 
this has been done by drawing the graph of o) for arbitrary (a)'s and altering them 
or adding new ones until the maxima and minima of the curve were all sufficiently 
small. Figs. 1 and 2 are graphs of w, representing two approximation processes 
requiring equal amounts of arithmetical labour. In fig. 1 the (a)'s are distributed 
over a wide range. "^ In fig. 2 all seven (a)'s are made equal to Xl^ The curves show 
that these distributed (a)'s reduce the amplitudes corresponding to a wide range 
of X^ to less than one-tenth of their original value. On the other hand the (a)'s 
concentrated at X^^ reduce the amplitudes in the neighbourhood of X^^ much more 
perfectly, but leave the (P)'s of graver period less affected. The allowable type 
of curve for co will depend on what is to be done with the integral of 2)^^ = when 
obtained. If its space-rates are required it is more important to abolish the modes 
of vibration having the larp:est values of X^ than it would be if volume inteerrals alone 
were needed. 

When the ratio of Xi^Xl^ is large, as in large tables, it is difficult to remove Pi by 
the processes indicated by figs. 1 and 2. For example if 3)^ = ^^/^x^-\-^^/^y^, 
and if the boundary is a square of ten co-ordinate differences side, on which 
^ vanishes, then the (P)'s having the lowest values of X^ will not differ greatly 

^ The advantage of distributing the (a)'s fairly uniformly was pointed out to me by Prof. 
A. E. H. Love. 

VOL. CCX. — -A. 2 T 
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from cos (-^Qmirx) . cos (Yomry), and for these l!^ ~ (m^+n^) j^-^tt^. So that the 
ratios of X/, X2^ &c., to X^^ when set m order of size will run 0*0247, 0*0617, 0*0987, 
0-1234, &a, up to nearly unity. Referring to %. 1 it is seen that the seven 
approximations would reduce the amplitudes of Pi, Pg, &c., in the ratios 4-0'48, 
+ 0*07, — 0*065, — 0*08, &c., the rest never exceeding -j^.^- and averaging about ±-^. 
Of course continued approximation would gradually reduce the amplitude of Pi, but 
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Curves illustrating the process of approximation. 



in cases like these it may be well to make a guess at the form of Pi, to estimate 
Xi^ as in Appendix, equation (33), to find an approximation to the amplitude Ai in 
4>i~'^u = SA^^Pa by the Fourier method (see Appendix, equations (29) and (22)), and 
so to remove the greater part of the first term of the series before the approximations 
are begun. This has been done in the problem of the dam, § 4. 

We have so far supposed X/, Xg^, ..., X/ unknown. If any one X/ of these be known, 
then making a = X/ will entirely remove P^ from the series. This process may 
sometimes be useful for removing the gravest modes of vibration Pi, Pg, &c. 

Since the value of <w is independent of the order in which its factors are multiplied 
together, it follows that the result of (j^t+i of a series of operations of the type 
</>m+i = <^m + ct^"'^ 2)'(/)^ depends on the initial guess <^i and on the values of the (ce)'s, 
but not on the order in which the (a)'s are taken. The application of this result 
to practice is slightly limited because the number of significant figures retained 
is necessarily limited. 

In carrying out an approximation-process with a set of (a)'s designed to make 
CO small, it is frequently only at the last stage that the predicted improvement 



BY FINITE DIFFERENCES OF PHYSICAL PROBLEMS, ETC. 323 

appears in the table. The intermediate stages may be wildly irregular if the values 
of a decided upon are used in ascending order of magnitude. If, on the contrary, 
the descending order is pursued the table tends to improve more regularly. 

§8 '2 '2. The Equation ^'(^ = p, where p is a completely known Fimction of 
Position — This can be treated in the same way as ^'</) = 0, except that the 
approximations must now be of the form 

for, since ®'(^m = p, this may be written 

which is the same as (3) of § 3'2'L 

§ 3 '2 '3. The equation 

(3)'-X^)P = (1), 

together with s homogeneous boundary equations such as 

where /3i ... ^^ are the boundary values, xjji ...xfj^ the body values and the (/^s are 
given numbers. P is now written in place of (^, because by (4) of § 3'2'1 P is defined 
as satisfying equation (l) of this section. We will suppose that both P and \^ have 
to be determined. As equations (2) contain no terms independent of the (i//)'s, it 
follows that when these expressions for the (y8)'s are substituted in the body-equations 
the latter become homogeneous, and are only consistent for the particular values of 
X^ which we have already denoted by Xi^, Xg^? • • • ^n- Further, on account of this 
homogeneity, any multiple of an integral satisfies the correct boundary conditions. 
Such a possibility does not arise with 2)'<^ = 0, for as 2)' contains no adjustable 
constant such as X^ it has no integral save (/> = 0, unless the boundary equations 
contain a term fp independent of the (t//)'s. And if they do contain such a term, any 
multiple of an integral fails to satisfy them. Thus, in the theory of membranes, if we 
put T) = 373^^4-3782/^, and (f) for the small displacement of the membrane from a 
fixed plane, then in the case of the membrane at rest 3)'<^ = 0, and <^ is commonly 
given at the boundary, its values there being /lo? ^20? • • • /so- 

On the other hand, for the vibrations of the membrane (2)^ + X^) <^ = the boundary 
condition is commonly ^ = 0, corresponding to the vanishing of all the (/)'s. 

To return to the general form of 3)^ : let P;^ be the integral desired ; (P/c)i, (Pk)2, • • • 
(P;^)^+i the initial guess and successive approximations to P;t- The fact that any 
multiple of an integral is itself an integral allows us to put for the body-points of 

(P/c/m+i = y(P;c)m~~^ 2j {^ k)m » . . f . r , , (3), 

2 T 2 
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which is more general than (2) of | 3*2*1 by reason of the factor y. The boundary 
values of (Pk)m+i ^^^ filled in so as to satisfy the boundary equations. Now (P;^)^ 
can be imagined as expanded in the unknown series 



\-*- kjm 



SljjX j 



(4). 



On the diagram of m and X^ (see fig. 3) equation (3) means that the straight line 




Fig. 3. 



representing a single process of approximation may now be drawn in any way instead 
of having to pass through m = 1^ X^ = 0. 

By suitably choosing lines all of which pass through ca = 1, X^ = X/, we can reduce 
the amplitudes of every P in the series except that of P/,, which is left unaltered. 
See, for example, fig. 3, where X/ has been given the particular value I^X^l To choose 
these lines we must know X/ at any rate approximately. For the first step an 
approximation to X/ is therefore calculated from (P;,)i in the way described in the 
Appendix equation (33). Denote it by (X/)i. For the succeeding steps (X/)2, &c., 
are calculated similarly from (P/,)2, &c. By Appendix equation (33) the errors in F 
are reduced more rapidly than those of P. The success of the method will depend 
on the original guess (Pk)u when expanded as tAjPj, being free from (P)'s having X^ 
nearly equal to X/. 

§ 3 '2 '4. Error in the Integral Due to Incomplete Approximation.- — A general guide 
here is the approximation process itself If, for example, this has been such as to 
diminish the amplitudes of all the Fs to less than yo ^f their former values, and if, 
for all that, ^ has not changed by nine times the permissible error, we may conclude 
that the process has been carried far enough.^ 

"^ Apil, 1910. — This is probable* but not certain. Thus if 99 sin x + 100 sin (3x) become^ 
9 sin aj + 10 sin {$£) the value at a: == Jtt does not change but it is not therefore ^ero, 





Table. 


<^1 


y'h i^/<^i 


h 


^ 'h ^ 'h 


h 


&c. 



BY FINITE DIFFERENCES OF PHYSICAL PROBLEMS, ETC. 325 

Again, from the size and distribution of 2)^<^ combined with a knowledge of various 
integrals of 2)^(/> = /{^^y), ^ rough estimate of the errors in ^ can frequently be made. 
Again for certain equations'^ the method of contour integration applied to a circle 
affords a check on the value of (j) at its centre. This method is very rapid, and it is 
particularly advantageous when applied to a circle enclosing many body values, for 
then a repetition of the approximation process would be correspondingly tedious. 

§ 3 '2 '5. Routine of Approximation. — Time and Cost. — To anyone setting out on a 
problem I offer the following experience as a guide in forming estimates :— ^It was 
found convenient to enter certain stages on a table with 
large squares, each divided into compartments. Thus for 

'^~l + ^ ^ 2^ 2 "^ ^04 = 0? C)ne 01 the squares is shown m 

the annexed table. All the quantities in it refer to the 

central point of the square. The intermediate stages are 

done on rough paper and thrown away. So far I have 

n 
paid piece rates for the operation 'S/+'S/ of about t-q 

pence per co-ordinate point, n being the number of digits. The chief trouble to the 
computers has been the intermixture of plus and minus signs. As to the rate of 
working, one of the quickest boys averaged 2,000 operations 'S/-i-'S/ per week, for 
numbers of three digits, those done wrong being discounted. 

§3*3. Relative Merits of Simultaneous Equations and of Sueeessive Approxi- 
mation. — The method of simultaneous equations may be applied to differential 
equations of any order and degree. It gives results which are exact for finite 
differences. It is necessary in discussions as to the existence and properties of the 
integrals of difference equations. But for actually calculating the integrals the 
labour becomes very great as the number of unknowns increases, and is of a sort 
which a clerk will not easily do. Large numbers of digits have to be dealt with, and 
a single mistake generally throws the result altogether out. 

The successive approximation methods of § 3*2 have only been applied to a limited 
class of linear equations. The results are not exact even for finite differences. But 
the bulk of the work can be done by clerks who need not understand algebra or 
calculus. Small and infrequent mistakes, or taking only a small number of digits, do 
not prevent one arriving at a fairly correct result. Nevertheless, it has been found 
best to have everything worked in duplicate. 

The method of successive approximation to the surface z = (^^ =y(^?2/) reminds 
one of the manufacture of plane metallic surfaces. The initial form of the surface is 
arbitrary in both cases. The essential things in both cases are a method of testing 
the work at any stage, a tool with which to alter the surface and judgment in using 
it. Methods of testing the arithmetic have been described in § 3 '2 '4 above. Our 

^ These include V2(/) == 0, Vi^t/, == Q. See a paper by BoGGio, * Jahrb, Fortschritte Math.,' 1900, p. 740, 
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tools (f)^^i = (j^m—^mT^^^^^m ^^^ tbo Founer method of removing principal modes of 
vibration may be compared to automatic grinding machines. The use of the hand- 
scraping tool corresponds to the adjustment of the numbers according to the 
judgment of the operator. This is always permissible. 

§ 3*4. The treatment of infinity may be illustrated by the following hypothetical 
example : — Suppose we were to set out to determine the field due to a rotating mass 
of gravitating fluid of known shape, the gravitational potential would have to vanish 
at infinity. But to be determinable by the methods of § 3 the potential would have 
to be ffiven over some boundary which could be included in the sheet of paper upon 
which its values were written. We might, for example, assume the potential at the 
edge of the paper equal to that due to the given mass of the liquid concentrated at 
its given centre of gravity, and find the figure of equilibrium on this hypothesis. 
(See treatment of base of dam in § 4.) 



§ 4. The problem of the determination of the stresses in a masonry dam has been 
chosen as a final example for two reasons : (1) Its practical interest arising out of the 
great expense and productiveness of these structures and the destruction of life and 
property should they burst ; (2) Its simplicity in that we have to determine a single 
quantity )( as a function of two co-ordinates only. 

The methods developed are suitable for finding the stresses, not in dams only, but 
in a prism or cylinder of any section, acted upon by any distribution of surface stress 
which is normal to the axis of the prism, when the surface stress is given, provided 
that the shifts parallel to the axis are zero. § 4*1 is devoted to theorems relating to 
any shape of contour. In § 4 "2 the special contour of the dam is introduced. 

The discussion that follows is founded upon two papers in the Drapers' Company 
Kesearch Memoirs (Dulau and Co.) : 

(1) '' On some Disregarded Points in the Stability of Masonry Dams," by L. W. Atoherley and 

Karl Pearson, F.E.S., 1904. 

(2) *'An Experimental Study of the Stresses in Masonry Dams," by Karl Pearson, F.R.S., and 

A. F. Campbell Pollard, assisted by C. W. Wheen and L. F. Richardson, 1907. 

I also owe some ideas, e.g,, the use of equipollent loads in the base, to suggestions 
thrown out by Prof. Pearson in the course of conversation. In these papers the dam 
is regarded as a prismatic-shaped body of indefinite length, so that the problem may 
be discussed in terms of two co-ordinates x and z lying in the vertical cross-section 
of the prism. The sluices are ignored — a serious omission. 

§ 4-l'l. Prof Pearson and his collaborators lay much emphasis upon our ignorance 
as to the real conditions at the base of the dam and the consequent inapplicability of 
analysis based upon special assumptions at the base, such as linear, parabolic, or 
quartic distribution of shear. Admitting this ignorance, there seems to me one 
assumption more reasonable than any of the others, namely, that the dam may be 
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regarded as a ridge upon the surface of a very large slab of rock which, with the 
dam, forms one homogeneous isotropic elastic system ; together with the further 
assumption that this substratum sufiered from no internal stress before the dam was 
built, except, of course, that due to its own weight. The consideration of the results 
that follow from these assumptions will be continued in § 4*1*9 below. 

§ 4*1*2. Cofiventions as to Symbols. — The Z axis vertically downwards. 

The X axis horizontally directed from the water towards the tail. 

0, the angle which a line makes with the Z axis, to be reckoned positive when the 
rotation is from the Z axis to the X axis through 270"^ and so onwards. 

The stresses xx and zz positive when they pull adjacent portions of the material 

together, and therefore the shear xz positive when the lower portion is pulling the 
upper towards the tail. 

q and .s distances along the outwardly directed normal to the masonry and rock 
and along the boundary drawn to the right of it. 

N and T the normal and tangential stresses on the surface related to q and s in the 

same way that zz and xz are related to z and x, 

I and n the cosines of the angles which the outwardly drawn normal makes with 
the X and Z axes respectively. 

§ 4*1*3. Position of the Origin, — Pearson and Pollard (p. 37) take this at the join 
of the front and the top. It will be more convenient in what follows to take it 
vertically below this point at the level of the surface of our hypothetical slab of bed 
rock. This is done throughout. The water surface stands at z = —p, and the 
pressure due to it is accordingly gp^ {-^p + ^), where p^ is the density of the water. 

§ 4*1*4. Specific Constants of the Masonry and Rock, — Following Pearson and 
Pollard I take the density p = 2*25 times that of water (p. 29) and Poisson's ratio 
Tj as |- (p. 33). 

§4*1*5. Units. — Distances are reckoned in metres. Forces in metric tons, each 
equal to the weight of a cubic metre of water. Consequently gp^ the weight of unit 
volume of water is equal to unity. And gp = 2*25. 

§ 4*1*6. Equations to he Solved. — ^The stresses may be expressed in terms of a single 
scalar y^ which satisfies 

^■*^ = (l=+^3S%' + l)x = « « 

S=^-S.p.. S = -^, S=-^. . (2),(3),(4). 



The surface equations then become 



av 3'- 



— p — z = — ^ ^— e 

dz^ dx dz 



_e(p + ^)= ---X+/^_grp2;U . ...... (5) 



dx dz \ dx 
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under the water. On the air surface the left-hand sides of these are zero. Here e is 
the tangent of the angle which the surface, drawn to the right, makes with the 
z axis. 

We shall require the surface stresses in terms of the space-rates of ^ along the arc 
and normal The theory of stress functions is fully discussed by Mr. J. H. Michell 
in^Proc. London Math. Soc./ Vol. XXL; on p. 110 we find the required trans- 
formation, which expressed in our notation reads as follows : — 

These equations refer to a weightless solid, hence the suffix. To complete them 

we must add the stress zz = —gpz^ xx = = xz after transforming it to components 
about the normal and arc. This gives 

§4*1 7. The Size of the Dam. — Suppose we have determined ;>( for a dam of a 
particular shape filled to a certain fraction of its height with water. Let us say 

Xi=/(^,^). 

We wish to find the stresses in a dam b times as big every way and containing 
b times the height of water. Try ^2 =/(^/^5 ^/^)- Then d\2l^x^, 9^2/9^^? and 
3^2/9^ 3^ would all be 1/6^ times their former values at corresponding points of the 
surface, but the constant terms in the surface conditions, due to the water pressure 
and weight of masonry, are now b times their former values at corresponding points. 
Consequently b^f(x/b^ zjb) is the form of ^ appropriate to a dam b times the size of 
the one for which f{x^ y) was determined, and the stresses in the former will be 
everywhere b times as great. For convenience in calculating, we will suppose that 
the co-ordinate difference, hx = Sz, is equal to the unit of length, unless otherwise 
stated. The result can afterwards be applied to a dam of any size. 

§ 4*1 '8. Simple Transformation Concerning the Density p of the Masonry. — If the 

reservoir be empty, and v^ is the integral for density pi, then ™ ^, is the integral for 

Pi 

density />, for it still satisfies the body equation V^ ( -^ Xe ) = ^ and also the surface 

equations (6) and (7), since both N and T vanish, so that the calculation of one 
integral suffices for all sizes and densities, if it be multiplied by the proper constants. 
Now when the reservoir is full we may take account of any density by means of 
two independent integrals. Let us calculate x^ f^^' reservoir empty and density pi 
and X f for reservoir full and density p2^ then the proper stress -function for reservoir 



BY FINITE DIFFERENCES OF PHYSICAL PROBLEMS, ETC. 329 

full and density p will be v^— — y^ + PX" ^^^' ^^ account of linearity, the solutions of 

Pi 



^^2 



1 a 



^ — r^ -^ — fi^gpz = Ni, = N2, = N3 in turn 



3^^ E, dq 

when added tos^ether satisfy ~~^~- ~~ -^ —n^gpz = N1 + N2 + N3. Here Ni is the 

normal stress due to the water pressure of a full reservoir, and Ng and N3 are the 
same when the reservoir is empty, that is zero. The above is the method of Peaeson 
and Pollard (p. 28) translated into stress-function symbolism. It shows us how to 
find the stresses for a full reservoir sustained by masonry of any density by 
calculating in detail two cases only. 

In this paper the case of full reservoir and density 2*25 is the only one treated. 

§ 4*1 '9. Integration of Surface Equations, — We see from equations (6) and (7) that 
if the shape of the boundary and the stresses upon it are given, then starting at a 

point 5o, and assuming initial values of )(, ^x/^S"' '^^^ ^x/^^> ^^ ^^^ ^^^ ^ double 
row of values of ;>( all round the boundary by straightforward integration. The 
initial values of ^, Sx/^?? ^^^ ^x/^^ ^^^ ^^^ significant, for they depend only on the 
arbitrary function Ax + B2; + C, which may be added to any distribution of x without 
affecting the stresses. The only outstanding uncertainty is at a sharp corner where 
1/R becomes infinite. Prof. Peahson has shown that at a sharp re-entrant corner 
the stress may become infinite (' History of Elasticity,' vol. ii., § 1711). In view of 
this it appeared possible that the boundary strip would be indeterminate at the 
corners, and the following investigation was made to settle the question. As it is a 
question of infinity, the finite forces N (i6', Tds^ and —gpzdsdq in the immediate 
neighbourhood of the corner may be neglected, leaving simply 

..... (8), (9). 

Now let us suppose that the trace of the boundary on the y plane is a circle at the 
corner, with the intention of making the radius indefinitely small after integration, 
on integrating equations (8) and (9) with R constant it is found, after some work, 
that 

-^ = Ce""K sinf — +D] 4-H (10) 

is the complete integral, C, D, and H being arbitrary constants. 

Now suppose that s/R, the angle turned through, increased from to a. Then it 
can be shown from (10) that 



B\ 1 ax 


3\ 


1 3x 


ds' E dq ' 


3^ dq 


R 8s ■ 



COS a+ (--A.) sin Gt, (3^) =( 5^) cos a— (-7/^) sin a . (11), (12). 



\ds /a \ds /o \^9^/^ \dq/a \d<l/o \ 3.9/0 

VOL. CCX. — A. 2 U 



330 



MR. L. F. RICHARDSON: APPROXIMATE ARITHMETICAL SOLUTION 



A 



I) 



It is remarkable that the changes in dyjos and 3x/3^ depend on the angle turned 
through, but not upon the sharpness of the corner, as is evidenced by the absence of 
li from (11) and (12). These equations are the conditions that the tangent planes of 
the surface x =/(^? v) j^^t before and just after the corner shall be parallel to one 
another. Thus the boundary strip of ^ can be integrated independently of the body 
equation, ^ ^^ = 0* 

§ 4*1*10, The following theorem about the total changes in dy^jdx, 9x/9^j and x over 

any length of boundary affords a useful 
check on the integration of equations (6) 
and (7). Let AB be any part of the 
boundary. Cut off in ima^fination a portion 
of the solid by horizontal and vertical lines 
through A and B meeting in D. It will be 
assumed at first that AB does not intersect 
AD, DB, except at A and B. Now apply 

to AD a normal stress zz^ and to BD a normal stress xx^, arranged. in amount and 
distribution so as to balance the forces acting on the boundary AB, together with the 
weight of the portion ABD. Let the weight be gpv acting in a line distant I from D, 
and let the stresses on the real surface AB be equivalent to a force X, Z acting 
through D, together with a couple G, The forces X, Z are to be reckoned positive 
when they are directed out of the solid. As the body ABD is in equilibrium x> 
dx/Sx, 3x/^^ ^^^ single valued, and if there are no infinite stresses they will be 
continuous. The changes along AB can, therefore, be obtained by integrating along 




rA 



AD and DB. To balance the horizontal and vertical forces —zZodx+Z+vc/pv = 

Jd 



and 



•D 



B 



lAytAJ. 







dz+'X = 0, where f is 4-1 if AB is above AD and —• 1 in the reverse 



condition. 



Now, substituting 






d' 



dixf 



X- 



9P 



z 



and 






-^^ , there results 

8r 



J.X. 

x + Z-^pgpv and 

D 



D 

dx/(^z = X, Also,, since the forces on AD, DB are 

B 



^^ 9x — 

D OX 

purely normal, d^l^z is constant along DA and 3x/^^ ^® constant along BD. There- 
fore, remembering that by (11) and (12) dxj^x, dx/dz are continuous at the corners, 
we see that 



B 



JL = —Z+gp{pv—Z()a) and 

A ox 



"9X 

aBz 






( 1 d K ( 14 ), 



where ci is the length DA. These are the required total changes in dx/^oc and 3x/3^. 
Lastly, to balance the moments about D, we must have 



•A 
D 



^A^Q 



iX — Xq) CtX~T' I XXq \Z' 

Jb 



Zo)dz+G-i-fJLlgpv = 0, 
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OO X 



where Xq, Zq are the co-ordinates of D, and where /x is + 1 if AB is to the left of BD. 
and — 1 if they are reversed. Integrate by parts, then 



D 



\^X XqJ 



/VvvA KA/*Xy 



A/'Z'i 






+ 



\% '~~ /Cqj XXq (a/%"™' I XXq KiX'/Cf) 






D 



a 



~\-G + filgpv = 0. 



Substitute the values of the stresses and 



B 



A 



X 



ax /A 



a' 



Hi/B+^H'^i-"'^^"^' 



where h is the length BD. Or, since by (13^) and (14^) (-^1 ~ ('^) "~^' 



B 



A 



X 



ax J A 



hi-J^] +hX+gp(zo~-—[xlv 



dz/A 



2 



(150. 



a useful formula. 

When the boundary of the solid intersects AD, DB, as does the curly line in the 
fig. 5, then draw horizontal and vertical lines to form a 
figure enclosing the boundary. Apply the stresses to 
FH instead of AD and to LK instead of BD, leaving 
the rest of AFHJKLB unstressed. These stresses 
must balance the surface forces as before, and now the 
whole weight of the solid enclosed between AFJLB 
and the boundary. On account of the change in level 

from AD to KJ or BL the term —gpz in zz is different, 
and this term beinff integrated alons; FJ and LB has 
the effect of subtracting the weight of AFJLBD from 
that of the solid part. Consequently we may integrate 
along AD, DB and keep everything else as before, 

pB 

provided we replace -hvgpv by ±gp (z—z^dx^ where z follows the boundary curve. 




A 



pB 



Similarly in (15^ we must replace ydgpv by ±gp (z—z^ {x—x^ dx, 

J A 

§ 4*1 'll. 'In order to simplify the arithmetic, the surface of the dam of the form 
chosen has been represented by straight lines, either horizontal, vertical, or sloping at 
45 degrees. It is one of the peculiar advantages of successive approximation methods 
that a simple case like this comes in conveniently as the first stage of the solution for 
any rather different boundary. 

Table of surface equations when the surfaces are horizontal^ vertical^ or at 45 degrees 
(fig. 6), deduced from equations above : — 

2 u 2 
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^'■X 



X 
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'^X^-^z-zsk-'^o 
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o 



~' - • - — =«• o 
'3x^ 



#^ /'^ 






OX^'l. 
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€-0 



C-0 






^=.->/i 



Zxd'^ 



o 




£ «« ± CO 



■HXac^ 



•^^ ^ ORIGIN 



BOUNDARY CONDITIONS 




^^""^ 



WMEN the: surfaces are: horizontal , vertical, or inclined at 4-5" 

Fig. 6. 

The rather complicated conditions on the sloping surfaces require special reduction 
to finite differences. Suppose the flank runs through a set of values of ^ ^^ ^^ 
Table VI., then at the point "5, 

d\fdxdz may be taken as i-{(Ol) + (11)- (ll)~- (01)}, 

d\/dz' as i-{(ir) + (01) + (11) + (01) -2(10) -2(00)}. 

Therefore the surface condition 

^'X + ^!X__==0 becomes = (Ol) -f (11)~(00)- (10) . . . 

Table VII. 



• — =- — p — ■ — 

dz^ dx dz 



(.13). 



Table VI. 
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In a similar way, at the point 0, *5 the second surface condition 







3^ 



)(._|.^X___2i^ becomes = (11) + (10)- (01)"(00)-2l^ . (14). 



dx^ dx dz 
The conditions on the front may be transformed in a similar way if necessary. 
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We proceed to show that the natural way of turning a corner in finite differences 
is consistent with (11) and (12). Take, for example, the corner in Table YIT. 
representing the tail rising vertically from the base. On the horizontal surface 

are taken along the line z = midway 



z 



'^2 r\2 

0, the equations — ^- = 0, tt-^- = 

3.x dx cz 



Table VIIT. 

Surface 



betw^een the two boundary layers. When we come to the corner we may either 

continue this process and so determine & = 5, or else we may take r;r 9 +Ma;M^ ;r-4~ — 

^x'^ dxoz 

at CC = +|- ; z = 0. 

Either process leads to ;>( = 3 at x = —^, z = —^, but one fixes h while the other 
leaves it arbitrary. Again we get x ~ -'^ ^^ ^ ~ '-^^ ^=—11"? either by taking 
d\/dz^-\-d\/dxdz = at x = 0, z =^ —^, or else by taking 6 = 5 and d\/dz^ == 0, 
and d\/dxdz = 0, as holding at the middle of the two columns on x =:-0. We 
may expect that when b is determined by the boundary equations a sharp corner is 
represented. When b is not so determined, a corner bevelled by the line x—z=^^. 
In the absence of the foregoing analysis 
these ways of turning the corner would 

have seemed tempting but risky. But we Surface ^ / 

shall see that they really do correspond to / 

the analytical way. For, on the horizontal 

bed, we have {^x/^^)^ ~ 1? (^x/^?)o — ""2, 
and the angle turned through is +-J7r. 
Therefore, from (11) and (12) dx/^s and x^— 
^xl^9. ^^ ^'^^ ^'^^ above the corner should 
have the values —2 and —1 respectively. 
And these are exactly what we find in the 
finite difierence table for the first differences. 

Again, suppose the tail rises at 45 degrees to the base. Starting along the 
horizontal the value ^=:3at,x= —l^z= — |- may be obtained by assuming (14) to 
hold at the origin, that is to say, at the corner. 

The other numbers on the sloping surface were obtained by using (13) and (14) 
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alternately, gpz being left out. Now we have 



ds 







''(l)."-^'""^- 



SO 



that equations (11) and (12) give us —0*707 and —2*121, respectively, for 
^ and -^ on the slope just above the corner. And these are identical with the 

first differences of the table when Ss and 8q are given their proper values of ^2 and 
1/^2, respectively. 

§ 4*1 '12. A Possible Expermiental Solution. — Now that we have shown how to 
integrate the boundary conditions, the analogy with thin plates will help us. For it 
is known (Love's ' Elasticity,' 1906, §313) that for a thin weightless plate, originally 
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plane and unacted upon by surface forces except at the edges, the displacement w 
normal to the plate satisfies the equation = ^r—r + 2 r, o^ o + -^r-r ^ Vj V, the ed^e 

conditions being such as to make the displacement purely normal to the undeformed 
plane and x and y being co-ordinates in this plane. So that if a sheet of thin steel 
" ferrotype" is taken and fixed at the edges in such a way as to make a double layer 
of values of w there equal to the integral of the boundary conditions of the dam, then 
the displacement w elsewhere would also be equal to the integral ^ ^f "the complete 
system of equations for the dam. The principal stresses in the dam are then equal 
to the principal curvatures of the plate. For a method of measuring them, see two 
letters to * Engineering/ October 25 and November 1, 1907. 

Whenever one has to solve the equation Vi\ ~ ^? ^^^^ form assumed by a piece of 
postcard bent in the fingers will be worth considering. As an accurate experimental 
method this would have the advantage over those of PEARSOisr and of Wilson and 
Gore, that displacements are applied instead of forces, and that the resulting 
displacement to be measured may be larger than theirs were. 

§4*1 'IS. To return to the Cofiditions at the Base of the Dam, — Prof Pearson 
pointed out to me that the effect of the stresses in the base of the dam, at a distance 
where the base of the dam subtends but a small angle, will be the same as that of any 
other statically equivalent system over the base, in particular to the force-at-a-point 
found by compounding the pressure of the water acting through the centre of pressure 
with the weight of the dam acting through its centre of gravity. We may therefore, 
in imagination, remove the dam, leaving a horizontal plane with this force acting at a 
certain point, the pressure of the water in front and no pressure behind. Suppose we 
had the distribution of ^ corresponding to this system of surface forces. Then the 
corresponding double set of values of ;)( at a considerable distance from the dam 
will differ exceedingly little from the true values, and if we keep them fixed and write 
in the upper surface values on the dam we may adjust the numbers inside by successive 
approximation to satisfy V\ = 0, and the result will then be exceedingly close to the 
true integral in the neighbourhood of the dam. This is what has been done. 

From the linearity of V^ and of the stress equations, it follows that if the stresses 
corresponding to a number of solutions of V\ = 0, when added, give the true stress, 
then the stresses derived from the sum of all these solutions will also be correct. 
With the aim of providing distributions of ^ in the bedrock which shall enable 
engineers to solve '' Pearson's Dam Problem'' for any shape of dam boundary, I have 
considered the actual distribution of ^ in the bedrock as made up of the three 
following parts, which must be added in the proper proportions as described below : — 

(i.) The term ^i — '--rs9P^^ = ~i^^' This gives Sxx = —gpz^ and so combines with 
the stress due to the weight of the bedrock to make s^ = 0. 

(ii.) In a weightless bedrock. The stress function r/zg, due to a point force in any 
direction at the orio-in. 
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(iii.) The stress function i/zg in a weightless medium, when on the surface z = 0, 
xz = = zz from x = to x = + oo, and zz = —p from x = Qtox=— co. 

The stretch s^ is made to vanish at .t = 4- oo by making Zxx = zz = —p in that region. 
At X = ±00 the stresses due to (ii.) will vanish, and therefore so will also the corre- 
sponding s^. Consequently, if we add together any multiples of the (;)()'s specified in (i.), 
(ii.), and (iii.), we have s^ = Sbt x = + co when the rock has its proper density />. 
(iii.) is to be so placed that its origin is at the meet of the front and surface, (ii.) is 
to have the point of application of the forces placed where the resultant of the water 
pressure and weight of the dam cuts the plane z = 0, 

It will be necessary to consider distributions i/zg and xfj^ in some detail. We are 
indebted to Mr. J. H. Michell for both of them. (See ' Proc. Lond. Math. Soc./ 
1901.) 

First, i/zg. For a point force of unit magnitude acting on the straight boundary 
at the origin of the polar co-ordinates r, d\ The stress function ^ is given by 
\jj2 = —TT~Wd' sin 6\ and the force is in the positive direction of the line from which 6' 
is measured. The stresses are as follows : — 

"^^ -r'dff'^rdr" rir ' ^^ ~ Or^ ~ ^' "^^ " 3rVr 3^/ ~ V^^^)' i^^^ i^l> 

Next, i//3. The Stress Function due to the Lake, — Michell shows that if unit 
normal pressure be applied along a finite length of a straight boundary of an otherwise 
unlimited plate in which the sideways stretch vanishes, then the form of the stress 
function is 

,/,/ = (r^,^-r'V)/27r (22), 

where r, <^ and /, <^' are polar co-ordinates centred A and B, and AB is the initial 
line. He further proves that the axes of principal stress at any point P are the 
bisectors of the angle APB, and that if this angle is equal to a, the magnitudes of the 
principal stresses are 

— (a+ sin a)/7r along the internal bisector. ..... (23), 

— {a— sin a)/7r along the external bisector (24). 

We want the limit of \\f^ in the neighbourhood of one end B when the other end A 
is removed to an infinite distance. Let the origin of our co-ordinates x, z coincide 
with B. Then let A be at a very great distance x = —-t, r'(^' = -f 2;, r'z=zt^-x, 

1 

V''3^ = -r- {^"^<^— (^ + ^)^}- Now remove the infinite but stressless tz and we have 

^3^ = {r^4^—xz}/2Tr, <j6 being equal to d+^n in the notation of § 4*1 '2. On the upper 
surface d^/dx^ = d^fdr^. Therefore the traction zz is equal to -I- (^/27r. 
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It should, of course, be negative for a pressure. Therefore we will write 



i/;3= ^^{r^^xz} = 



l^!{x' + z')(d+'^]+xz 



(25), 



I have also arrived at this result by an independent method. 

We must next consider the stretch 5-^. For on the surface a great way in front or 
behind the dam one would expect s^ to vanish on account of the uniformity of the 
surface pressure. Now in Michell's solution for a finite loaded portion AB all the 
stresses vanish at an infinite distance, and therefore the stretches also. But when A 

moves off to infinity we find from (25), after differentiation, that [jlSj, = 
_l^i{2/9-4cos6^sin6>} 



8 



8 



Here s^ vanishes when 6 =■ — A-tt, that is, behind 



the tail. But not when = -^-^tt in front, under the water. 

§ 4*2. The detailed working for a partieular shape of eontour, 

§ 4'2'1. Tlie form ehosen for investigation is shown in fig. 7. It was the best 
representation of the Assuan dam, as drawn in Pearson's papers, which I was able 
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Fig. 7. Contour chosen for investigation. 

to obtain with so few as six co-ordinate difterences to the height, without inter- 
polation on the boundary. As a real structure it would be liable to crack at the 
points opposite the re-entrant angles on the flank, somewhat as the unfortunate 
Bouzey dam did, but that tendency will not affect the stresses lower down, with 
reference to which Pearson has given warning, and to which attention will here be 
directed. The height, _p, is 6 metres. The area of the cross-section above the line, 
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2: = 0, is ]3'25' square metres. And this, at a density of 2*25, gives a mass of 
29*81 metric tons above the rock surface per slice 1 metre thick. The centre of this 
mass was found by graphical construction to be at the point x = +1*335,^ z = —2*19. 
The total horizontal thrust of the water on the same slice is 18*0 tons. In dividing 
the forces acting on the structure into lake pressure and a force-at-a-point, after the 
manner of \ 4*1 '13, I have considered the lake to extend right up to the vertical face. 
This leaves a force of 0*125 ton acting upwards near the corner, due to the bevel 
diminishing the depth of water just there to be included in the force at a point. 

In finding the force at a point, the total horizontal thrust of the water of 18 tons 
acting through a point | of the height of the dam from its top, is compounded with 
the weight of 29*82 tons acting through the centre of mass and with the 0*125 ton 
acting upwards near the corner. The resultant of 34*74 tons strikes the base at an 
angle, the tangent of which is 1*657 (= 58"^ 54'), at a point the co-ordinates of which 
are x = +2*55, 2: = 0, If there had been a pronounced curve at the bottom of the 
front it would have been necessary to use a link polygon to find the load point of the 
horizontal 2; = 0. 

§ 4*2*2. Taking six co-ordinate differences in the height the surface conditions were 
integrated by equations (13), (14), &c. A difficulty arose at top. For the dam there 
being only one co-ordinate difference thick, that is, ;)( being expressed by three 
columns of numbers only, the four boundary conditions cannot necessarily be simul- 
taneously satisfied. It was avoided by making the total change in d^jdx over the 
base equal to the total downward force, and by making the total change in 9x/3^ over 
the base equal to the horizontal thrust of the water, and then integrating the boundary 
strips separately from in front and behind, and making y^ one-valued where they met 
at the top. 

On the lower boundary ^ was made equal to the sum of the following four parts :- — 

(i.) For the lake pressure, the values of 61//3, found directly from the co-ordinates 
X and z. They are given in Table IX. 

(ii.) For the force at a point 34*74\//2, found from the factors r and r sin 6^ measured 
on a large sheet of scale paper. 34*74i/;2 is also given in Table IX. They may be of 
use to anyone re-working the problem for a dam of a different shape. 

(iii.) A linear function added to the above two in order to make the sum agree 
with the values of ^ on the air surface behind the dam. It would have been simpler 
to have altered the surface values to make them agree with the base, but I did not 
think of that in time. This linear function is most conveniently specified by its first 
differences, which are d^jdx = +9*735, d^fdz =^ —15*35, and by the fact that it 
vanishes where the point-force cuts z ^ 0, These first differences w^ere calculated 
from the first space-rates of the expression — 34*747r"V^ sin ^ along and normal to 
z=i —58° 54', and were used to calculate the linear function. 

(iv.) Finally, —^z^ was calculated and added. 

* By calculation more exactly 1 * 340, 
VOL. OCX. — A, 2 X 
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The sum of these four 
functions of position is given 
on the lower boundary of 
Table X. Its lake and 
point-force components have 
been checked by plotting 
curves and recalculating any 
number which did not lie 
smoothly in the series. In 
this way an error of 1*0 
could be detected, and 
therefore probably does not 
exist. The multiplications 
were made with a 20-inch 
slide rule. Altogether I ex- 
pect the ^^ standard devia- 
tion" of the errors in the 
lower boundary values is 
0*15 or less. 

The initial body values 
assumed may theoretically 
be arbitrary, but of course 
the approximation process 
need not be so long if they 
are near to the correct inte- 
gral. In this case the initial 
body values were made equal 
to the sum of the four 
functions 6i/;3, 34*74i/;2, the 
linear function of (iii. ) above, 
and —^z^ in the region below 
z = 0. The numbers were 
read from graphs in a very 
rough way. Above z = 
the initial body values were 
simply guesses. 

The preliminary stages of 
approximation were full of 
experimental processes which 
did not completely succeed 
^nd are not worth describing, 
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especially as the arithmetic was somewhat faulty. Let us start with the table as left 
by them and call it ^i. 

Next the processes recommended in § 3 '2*1 were adhered to. The table being large 
it was thought best to remove Pi and Pg from ^i— ^,, = tA^Vf, by making guesses at 
the form of Pi and Pg, calculating their coefficients Ai and A2 by the Fourier process 
and then subtracting A1P1 + A2P2 from ^1 — ^^. Of course P;^ is defined so that 
SIP/ = I5 where S denotes a summation over the body-points, whereas the guess is 
an approximation to B/tP/^, where B is an unknown constant. For 2)' = ^ ^, and for 
zero boundary values of the (P)'s the general equations of the appendix give 

1=1, V = V/T = S [ ^ ^ (B.P,)]VS [B,P J. 
Now ^ '4>, = tK'^.'P,. Therefore S [P, . ^ *<^ J = X,'A,. Therefore 
A,P, = P, . S [P, . ^ *<^ J X S [B,P J/S [ ^ ^ (B,P,)]^ = (B,P,) . S [B,P, . ^ Vi]/S [ ^ ' {B.P,)]^ 

In the right-hand side P;;. only appears in the combination H^JP^c, so when this is 
known A/,P/, can be calculated at each point. A rough check on the accuracy of the 
guesses at BiPj and B2P2 may be obtained by comparing the values of X/ and Xg^ 
obtained from them, namely, 0*18 and 2 '3, with the corresponding quantities in 
infinitesimals for simpler geometrical figures of the same area. The area of this 
smaller dam table is 113 square units. And for a rectangle f as broad as long the 
formula of § 3*1 gives X/ == 0*15, while for a circular plate clamped at the edges the 
principal vibration with a single circular node has X^ = 1 '2 (Rayleigh, ^ Theory of 

Sound/ S 221a). It is possible that what 
has been called X2^ above was really Xs^. 
Finally the table was given four approxi- 

^ ^ mations of the type X»+i ~ X»"~^~^^^X» 

where a was 10, 30, 50, and 64 in turn. 

^ ^ These numbers were chosen because they 

V ^ gave a good *' curve of ratio of final to 

^Ss^ initial amplitude. '' The result is shown in 

S^ ^\^ Table X. From it the stresses can easily 

-"-^■^ -'^— — be calculated. The values of p'^x which 



X 



X a 



should vanish are given there also. 

I 4 '2 '2*1. Interaction of Body and Sur- 
face Equations. Bevelling of Re-entrant 
^ ^ X X Angle at Front—It has been shown in 

§ 4*1*9 above that the given stresses on 
the upper surface determine a double layer of values of x covering this surface. 
At almost all points of the double boundary layer it is impossible to evaluate ^ \ 
because one or more of the values of x involved is lacking, so that x over the 
surface must be determined by the surface condition only. The only exceptions are 
near a re-entrant angle. For example, at the points where x = ^ ^i' ^ 1^ ^^^ 
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Table XI. sufficient surrounding numbers exist to determine ^\ completely. 
We might therefore determine a and b at these points by the body equation. 
But then we should be left with a boundary strip indicating a considerable 
amount of bevelling. Or we may go to the other extreme, and determine ^ in 
the chequer marked c by assuming the surface equation ^\/^x^z to hold at the 
angle. I take it the boundary strip would then represent a perfectly sharp right 
angle. On these grounds I suppose that the case actually studied, in which a and b 
are determined by the surface equations and c by the body equation, represents an 
angle with the slight amount of bevelling caused by joining the points x = 0, z = — |-, 
and ic = — |-, 2; = 0. This was not realized till about the 10th approximation, when 
Mr. BoRCHAED pointed out some inconsistency. It was corrected along with other 
errors in the boundary conditions, and all final numbers refer to the bevelled angle as 
just stated. 

§ 4 '2 '2 '2. Error due to a Point-Force having been substituted for the Actual 
Distribution,' — ^We see from Table X. that the actual distribution of stress differs from 
the assumed point-force in being spread out well over the base. We may form some 
estimate of the order of the error involved by comparing the stresses over the lower 
boundary due to a vertical point-force with the same due to a statically equivalent 
pressure spread uniformly over the base. Miohell's stress function ^Q^{r^^—r^^^^)/2iT 
enables us to do this. (See equations (22), (23), (24).) For when the point P on the 
lower boundary is vertically below the centre of the stressed surface, the total force is 
proportional to the stressed area, which is equal to 2r sin ^a. So that if the total 
force is to be constantly unity as the area alters, the stress function must be equal to 

___ — _ . ___r_^ — XL, By (23) and (24) the principal stresses at P are 
/ir sm ifOL ZiiT 



rir 



,2 



. 1 — ---^^^^ + terms in a^, &c. I , ~ 

4.0! J TTT 



.2 



—- + terms in a^, &c. 
3 ! 



In the limit when the force is at a point we have a = 0, and these reduce to -~2/r7r 
and zero. In the case of a point on the lower boundary, eight units below the level of 
the rock surface, the base of the dam subtends an angle of about 4|-/8 radius. Then 
06^4 . 3 ! = 0*013 : that is to say, the errors in the stresses at the lower boundary, due 
to substituting a point-force for this statically equivalent pressure, spread uniformly 
over the whole width of the base, are about 2 per cent, of the greater of the two 
principal stresses at the lower boundary. This uniform spreading is not, of course, 
exactly what has happened in the approximation process, but it is sufficiently similar 
for the question at issue (see Table X. ). The correct lake stresses and the correct 
upper boundary will tend to swamp this 2 per cent., which is therefore quite negligible. 

§4 '2 -2 -3. Errors due to Incomplete Approximation, — The values of ^\ given in 
Table X. for the last approximation look as though they consisted chiefly of a principal 
mode of vibration which had a single nodal line sloping parallel to the flank of the 
dam across the middle of the table. For this mode of vibration X^ will be considerably 
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greater than for the gravest mode which has X/ = 0*15 ; it will probably be as great 
as unity. Then the errors in ^ must be roughly equal to ^ \, If we imagine the 
tabulated values of ^ ^x smoothed to eliminate all vibrations, except the one with the 
single nodal line already referred to, and if we then take second differences of the 
smoothed ^ \, these second differences will be the errors in the stresses. It is easy 
to see that they will be less than unity. 

§ 4 '2 '3. Halving the Differences,— -Some notion of the error due to having so few as 
4|- co-ordinate differences in the base of the dam may be formed from the error in ^^, 
found from the equation (^*— A;^)(^ = 0, for the gravest mode of a square of seven 
differences to the diagonal. This is shown in §3*1 to be 13 per cent. It was 
therefore thought desirable to halve the differences and reapproximate. Body values 
half-way between those of the smaller table were filled in by interpolation. 

Since the boundary formerly lay half-way between two sets of numbers, a set of 
interpolated values now lies directly upon it, and consequently the surface conditions 
resolve themselves into relations between the three outermost layers of numbers 
instead of the two outermost as before. 

§4*2'3'1. The Method of Approximation. — ^Affcer some preliminary experiments, 
the numbers just inside the boundary were corrected by taking ^^.^^ = X^'^^~^^\n, 
and after each approximation the numbers just outside the boundary were corrected 
so as to keep ^xl^9- ^qual to SxP?? calculated analytically. The following values of 
a were taken in turn: 10, 30, 50, 64, and again 2, 5, 20, 40, 40, 50, 60. These 
numbers were chosen, as all such have been, because they gave a suitable curve of 
reduction of amplitudes. 

§ 4 '2 '3 '2. An analytic integration of the boundary conditions was carried out, 
using equations (6), (7) together with (11), (12) at the corners, and starting from 
^ = 0, '^xl^^ ~ ^' behind the tail of the dam, as in the small table. The expressions 
deduced for x ^^^^ ^x/^^ ^^^ given in a schedule. Values of ^ calculated from them 
were set down in their places on the boundary, and the first difference across the 
boundary was made equal to ^x/^^- These numerical values may be seen in 
Table XII. They have been carefully checked at a number of points by the theorem 
about total changes given in §4*1 '9, and have been found to be correct to 0*01 or 
less. The exact values of x obtained from the boundary have been compared on 
z = with those found from the point-force and lake pressure in the bedrock- — 



X. 


calculated strictly along 
boundary. 


calculated from 
point-force assumption. 


-4-5 
-5-5 


270-30 
329-987 

■ 


270-20 
329-90 



a quite satisfactory agreement. 
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§ 4 '2 '3 '3. Various adjustments were necessary. First, because the coarse-difference 
table has on its upper boundary numbers calculated stepwise by finite differences, and 
these are not quite the same as the infinitesimal values. Consequently, the infini- 



U-83 



17-20 • ^6-91 #39-9S 



Table XII. 



10.36 17-20 



10-97 17-20 



\0'BZ l7-;20 



26'S3 
-0-3S 

Z7-00 

-0-06 

Z7'50 

+ 0-31- 



The stress function - x for a dam 6 metres 
high. Co-ordinate difference one-half. The 
boundary of the solid passes through the 
series of dots. The black line encloses the 
body region. The small numbers are 

- tV ^ *x- 



6 '06 10-40 



17-60 Zd-00 

-0-05, -t-O'OS 



3-57 5-36 
3-33 5*86 
"^-.'ZO 5-56 



10-60 18'fO Z5-19 

+ 0-'^Z -0-I6 +0-93 

(1-02 \d'5d 30-57 

-1-81 40"V-9 ~0'I6 

12-06 Zl'-fe 33-3^ 

-O'OI -O-'^-'i- -^I'OZ 



Z-25 3-00 



7-05 14-25 ;i4-37 36-67 

+ 1-23 -1-5.5 +0-77 tO-5l 



I -23 1 .27 



3-97 3-43 17-63 2H-23 ^0-05 

+ 0-BO ■♦■I-I4- +0-75' -i-0-X& +0-09 



0-30 0-37 



Z'l\ 6-25 IZ'ZZ Zl'^e 33-00 46-2,5 

i-0"^d -O'OI -0-37 tO-53 -1-77 +0'A5 



0*65 0-55 0-33 0-05 
• 000# O-OO* 0-00# 000« 0-00 



1-22 4-42 9-54 17-56 27-21 39-02 52-77 

+0-50 -O'^O -0-37 +0'4S- -0-S3 -0"t3 tO-OS 



0-47 0-57 



Z'Z5 Z-57 



5'\B 5-35 



3-58 10-^6 



0-97 3-67 5-22 I4'g0 2.3-37 33-53 46-17 60-25 75'5d 

"O-ZS ■fO'17 -0-S& -t-O-Sd -f-O-li -eO-l^ -0-97 -0-03 •tO'70 



39 -33 56-43 

4015 56-S^ 

# 
40-^4 56-54 

41-31 57-54 

42-55 55-73 

44-48 60-41 

47-13 62-58 

50-65 66-21 

35-17 70-49 
60-51 75-52 

67-7/ 5^1-90 99-54 

91 •55«105-64#l27-24#l4-7-33< 



0-65 0-55 (-67 ;3-54 7-75 13-52 2,1-03 30-44 41-57 5f43 65-75 84-25 100-53 |j7-54 

+0-/0 -0-67 +0-09 -0-3'V ~0-05 t03l -0-63 -0-U -0-55 -(•0-64- -0*19 +0'S5 -^O-SZ -0-32 

2-96 C5 -65 5-25 3-50 13-46 20-13 25-56 35-7) 50-35 63-55 75-07 33-70 110-23 127-57 

-f0-+6 "O'SZ -1-17 -0-52 +0-IS -0-44- +0*03 +0-69 +0-15 -O-a.5- -fO-37 fO-08 ^|.3.5 -tO'lS 

6-SL 5-35 10-55 14-95 20-65 2S-II 37-22 47-93 6007 73-53 55-16 103-31 120-62 135-25 

+0'I3 "O-IZ -1-5,1 +0'S7 -0-4-1 .+0'll -O-OI -O'CJ- -O'lS" +0-32 -0-36 -♦• 0-A.7 +0&Z +0M3 

12-47 14-77 15-13 22-95 23-23 37-3/ 46-96 55-14 70-67 54-35 33-12 114-93 131-73 149-56 

-0'I3 -0-S7 +0-53 i-O'lSr i-O-OZ -0-5.7 O'OO -to- 14- fO-70 +02,4- +0-5-3 -0-;i3 -0-33 -0'09 



17-24 19-62 



22-53 26-71 32-19 39-05 47-66 5771 69-26 52-05 36-01 110-35 l2£-99 143-35 

-0-09 -0-14- -0-13 -0-5.5. -004- +0-5.5, -0-50 -0-+8 - o-5,5 -0-37 - 0-43 +-0-IO 



136-57 156-60 



14613 166-73 



157-26 177-59 



168-43 183-28 



162-11 I51'22 



ZS-ZZ 32-33 36-43 42-65 49-9.5 59-09 69-fO ^hS7 3450 105-68 123-54 140-14 15713 175-71 



42-35 47-35 54-i;l 6I.-50 71-52 52*06 94-34 107-60 1,^2-26 13732 154-10 171-28 



tesiinally correct upper boundary of the fine-diflference table did not at first fit the 
lower part derived from the coarse-difference table. A correction to the coarse- 
difference table was therefore calculated, which, if added to it, would have made its 
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upper boundary infinLtesimally correct as regards ;)(, when the value of ;)( was found 
by interpolation^^ i^x/^^I appeared to be near enough already). And this correction 
was used to bring the lower boundary of the fine- difference table into harmony with its 
upper boundary. Its greatest value in the region of the fine-difference table was 1*5. 
Secondly, a principal mode of vibration, having maximum value 0*5 and X^, roughly 
1'5, was left in the top of the dam not completely removed by the approximation 
with a = 2. It was removed by guesswork satisfactorily. Thirdly, a number of 
slips in the arithmetic had to be corrected. The result of these processes is shown 
in Table XII. 

§4 '2 '3 '4. Errors Due to Incomijlete ApproximatAon.- — ^As a justification of the 
body values of Table XIL, I propose to consider the final distribution of ^/x 
belonging to them, and not the process by which they were obtained, for owing 
to errors and experiments this process was long and complicated. Let us imagine 
the distribution of ^I'^x gi^^i^ '^^ Table XII., to be expanded in the series 
tV ^ i*X — AiXi^Pi + A2X2^P2 + ... 4- ... + A^X/P^, where P/, is an integral of 
(i\ ^i"^— X/) P/, = 0, and has P;, — and ^P/,/'S (normal) = on the boundary. Then 
the error in x is APi + A2P2+ ... + A^P„. An inspection of ^i^x shows that it changes 
sign many times in a small area. Hence it obviously contains the principal modes of 
vibration for which X^ is large. The greatest value, namely, X/, is about 64, so for 
this the error in x will be -^-4 of ^6 ^ /x- The greatest value of yg ^ i*X ^^ ^^^^ ^^ ^^ 
1*81 5 corresponding to an error of 0*03 in x- If this had been distributed all over 
with alternate + and -- signs the corresponding error in the stresses would be 
(4x 0*03)/A^ = 16 X 0*03 = 0*48. This is appreciable, but not serious. The mean 
value of Ye ^ i^X' formed by squaring, adding, and taking the square root, would be 
much smaller than 1*81, and the error in. the stresses due to the presence of the 
higher modes of vibration P^, P^-i, l?n-2^ &c., will be correspondingly less than 0*48. 

Next, as to the gravest mode Pi : it is possible that Pi may be prominent in the 
error in x, and yet Xi^Pi obscured in ^/x ^J ^^^ presence of X/P,„ since Xi^ is a small 
fraction of X/, yVe can estimate it very roughly by means of the theorem that at 
the centre of a circle an arbitrary function f differs from an integral x of V^x = 0, 
which coincides with it, as to value and normal space rate, everywhere on the 
circumference, by the integral over the interior of the circle of the product 



•^ 167r 



^:yiog/ry+i-^(r"" 



= V^f X B, say 3 



where a is the radius of the circle. The integration has been effected with sufficient 
accuracy by drawing contours of B on tracing paper and laying it over Table XII., 
and then adding up the values of ^ /x situated between each pair of contours. 

■^ The internal distribution of this correction was calculated by contour integration applied to freehand 
graphs in a way which the author hopes to publish shortly. 
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The following results have been found : 
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* P is the excess at the centre of the tabulated number ~ x over an integral of ^ ] V = 0> which has 
the same value and radial space rate on the circumference, 
t Q is twice the previous column divided by square of radius. 

The corresponding errors in the stresses cannot be determined easily. As a very 
rough guide to them^ take the second differences coefficient at the centre of the circle 
measured along a diameter with co-ordinate difference equal to the radius, that is 

— (twice error at centre)/(radius)^. These 
differences are also tabulated and are seen to 
be small. 

§ 4 '2 '4. Peculiarities at the angles of 135 
degrees on the Upper Boundary. — There are 
three alternative ways of treating the relation 
between the body and boundary equations at 
the re-entrant angles. (See Table XIII.) 

(a) The number (lO) is determined so as 
to satisfy ^ iX — 0? ^^d does not enter into 
the boundary equations at all. (10) is derived 

from (01) by means of the known value of 
^x/^?' half-way between them. Also (11) is 

derived from (11) in a similar way. This is 
the procedure which has been adopted in calculating Table XII. 

(6) The number (10) is determined so as to satisfy ^i^x ~ 0? ^^^ ^^ ^^^^ enters 



10 



00 



10 



\ 



\ 



S 



01 



N 



11 






S 



s 



Body 
region 



\. 



12 X^ 

Boundary ^\ 



BY FINITE DIFFERENCES OF PHYSICAL PEOBLEMS, ETC. 347 

into the boundary equations, in that the average of (11) and (10) is derived from 

the average of (10) and (11) by means of the value of 3x/^? ^* ^^^ centre of the four 
points where these values are situated. 

(c) (11) is derived from (11) by means of the value of d^/^q half-way between 

them and (21) is similarly derived from (12). But (01) is derived from (10) by dx/^q 

at ^, ^ and (10) is also derived from (10) by d^/dq at 00, where in the latter q^is 

normal to the line joining 01 and 00. Thus (10) and (01) are both connected to (10), 
and are no longer independent. It is therefore impossible to satisfy the body equations 

completely. Large values of ^i^x ^^J ^^PP^a^i' at 10 and 01. In some work by this 

method, now rejected, the approximation was so arranged 

that these values of p" i^x "tended to become equal and Table XIV. 

opposite. The reason for choosing (a) was that the system (c) ^ 

gave large oscillations in the stresses near the corners, 

oscillations for which it was difficult to find any physical / 

reason ; (6) gave extraordinary oscillations in the shears on / 

the boundary, but (a) gave smooth curves of stress. X 

To confirm this choice the co-ordinate differences have been 
altered so that the boundary lies half-way between two 
tabular numbers. The ambiguities now disappear. Thus, at 

the hind toe it is found on stepping out the integral of boundary equations that we 
have zero values, as in Tables X. and XIV., whatever be the size of the co-ordinate 
difference. 

Then at the corner ~-^ = ~.^ _-A_ = - -~ . Therefore - r-^ h^r4r = 1. 

Now test Table XII. to see whether it gives the same result. 

At the hind toe |4 = +3*88, fi^fi^^^r = -379, and ^=|^ = 1*02, which is 

'die ^xm —3*79 

satisfactory. 

At the front toe a more elaborate check has been made by means of a new table 
with smaller co-ordinate differences, -^ = ^. This is Table XV. It was prepared as 
follows :— First, the number of digits in the values of x, near the corner, w^ere reduced 
by taking instead -^^ = -X"~''4 + 30.':C"~18;^, so that x and x^ have identical second- 
differences. Next, curves were plotted showing the variation of x' along lines parallel 
to the co-ordinate axes. Values of x^ were read from the curves at the new points 
required. These values provided the boundary numbers of Table XV., on the side of 
it which is in the middle of the dam. They also provided the initial body values, 
which are given in parentheses in Table XV. The boundary numbers on the 
masonry-water surface were, however, not obtained by interpolation from Table XII., 
but instead from the infinitesimal integral of the boundary conditions. Each of them 

2 Y 2 
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waB derived from the value of -y^ at the nearest point on the masonry- water surface 

by adding or subtracting the value of 7^ at that pointy multiplied by the distance to 

aq 

that point. The final body values of ^^ were arrived at by the process which has 
been compared to hand-scraping in § 3*3 above. They are justified by the very small 
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values of p^iV belonging to them. They are given to three places of decimals merely 
in order to get smooth curves of stress. The interpolated boundary values may well 
be ±0*02 in error, so that the third decimal has no real significance. The stresses 
deduced from this table are shown in two cases in figs. 8 and 9. It is seen that the 
sharp oscillations revealed in the stress by the fine difierences are passed over by the 
coarse difierences. The relation between the stresses near the front toe is, however, 
not consistent with the proportionality of the error to h^. This indicates that the 
assumptions which have been made as to boundary conditions have not been quite 
consistent with one another.* Nevertheless, the scheme adopted with h = ^ is to be 

"^ Or that the term in ¥ in the error cannot be neglected. 
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Fig. 8. Vertical stress zz on z = + J for a dam 6 units high. [Table (14) should read (15).] 
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Fig. 9. Shear a;^ on 5; = for a dam 6 units high. [Table (14) should read (15).] 
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preferred to the alternatives b and c, for with them the inconsistency was found to be 
greater. It would be interesting to rework Table XV., taking, on the masonry- 
water surface, numbers stepped out by finite difierences exactly as was done for the 
table with /^ = 1. It is not possible to adjust the boundary values of Table XV. so 

as to give by interpolation the infinitesimally correct values of ^^ and 7^- on the 

J. 

masonry- water boundary, for this would necessitate four independent adjustments in 
the three values nearest to either corner. 

At the protuberant angles of 135 degrees the method (a) adopted for the re-entrant 

angles will not apply. For according to this method number (10) is not determined 
by the boundary conditions (see Table XIII.). Neither is it determined by the body 
equation, for it is right outside the boundary. But it cannot be omitted, for it enters 
into the body equation at 1, 0. 

It has been derived from (10), as in the method (c). Fortunately a protuberant 
angle of 135 degrees only occurs at one point, and even if the method adopted is not 
quite correct, the stresses will only be affected in the immediate neighbourhood of the 
angle. This point is more than half-way up the flank. 

§ 4*3. Conclusion. — ^An attempt has been made to provide tables of the stress- 
function, giving by their second differences stresses sufficiently accurate to be of use 
to the practical designers of dams. The evidence that this has been done is : — 
(i.) The discussion of the boundary conditions in §4*1 '11, §4'2'3*1, and §4*2*4; 
(iL) The discussion as to the completeness of the apjjroximation in § 4*2*2*3 and 
§ 4*2*3*4 ; (iii.) The general agreement of the stresses derived from the three sizes of 
co-ordinate difference in figs. 8 and 9. In this statement about the accuracy it is 
assumed that the stress is taken from the table which has the smallest co-ordinate 
difference. That is, from Table XV. having h = ^ near the front toe ; from Table XII. 
having h = ^ elsewhere in the dam ; from Table X. having h = 1 deep down in 
the bedrock. 

An additional confirmation is the strong resemblance of the curves in figs. 8 and 9 
to those found by Messrs. Wilson and Goee"^ by stressing an india-rubber model 
with a rounded angle at the front toe. 

The arithmetic has been carried out by a number of people, principally Messrs. 
TiLLEY, H. BoRCHARD, W. Sheppard, C. H. Masters, and G. RoBmsoN, and I am 
grateful to them for the care they have bestowed upon it. 

§ 5. This work on the solution of physical problems by finite differences has been 
carried on at intervals during three years. My thanks are due to Mr. G. A. Schott 
for convincing me of the desirability of discovering new methods for solving physical 
problems; to Mr. A. Berry, Mr. G. F. C. Searle, and Prof Karl Pearson for 
their encouragement and advice in the early stages ; to Mr. A. Campbell, 
Mr. J. G. M. Garnett, and Mr. H. H. Jefecott for references and the loan of books ; 

* ^Institution of Civil Engineers/ February, 1908. 
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to Prof. A. E. H. Love for several vital improvements which were introduced on 
re-writing the paper, especially for the suggestion of the theory in the Appendix ; 
also to Prof. Larmor, Prof. Sampson, and Dr. B. T. Glazebrook for much valuable 
criticism. 

Appendix. — Properties of the Principal Modes of Vibration.^ 

In the theory of the approximation process we have assumed that — 

(i) An arbitrary function / can be expanded in the form / = %K.j^j,, where Vj^ is an 
integral of ±2)'Pa = X/P, and satisfies the same boundary condition as f does. 

(ii) That the X^ are real, all of one sign, and lie in a finite range. 

(iii) That a certain one-signed function of position I can be found such that 
SIP/,P^ = 0, SIP/ = 1, where S denotes a summation over the body points. 

It will now first be shown that if the system of difference equations satisfies certain 
conditions, then the simultaneous equations 

3)Vi = 0> 2)'t//3= 0, ... 1)'i/<„= of (1) 

are equivalent to 

91 = ^-0 9^-0 (2) 

where V is a one-signed homogeneous quadratic function of \fji..,xfjn. Next, it will be 
shown that the desired properties (i), (ii), (iii) can be deduced from the existence of V. 

No attempt will be made to determine whether the properties (i) and (ii) can hold 
under more general circumstances when V does not exist. J However, V exists for so 
wide a range of physical problems that it is well worth considering. 

Conditions of Existence ofY. — As equations (1) are to be equivalent to equations (2) 
it must be possible to find a set of numbers %, ig, ... 4 such that 

BY 
i0xlfi=:—-- for l=l,,.n (3). 

dxjji 
And then, since - is independent of the order of the difl^erentiations, 

^-pC^'^^'^O -= ^T-(^*^^'^^) • • ' • (4) 

o^k dxjji ' 

* The dynamical analogy which is the basis of this section was pointed out to me by Prof. A. E. H. 
Love. It has been introduced by Pockels in his book *Uber die Gleichung, /\H + ¥u = 0.' See also 
Eayleigh, * Sound/ vol. I, chap. IV. 

t The boundary values in equations (1) are supposed to have been expressed in terms of the body 
values by the boundary condition. 

X As to the property (iii), note that by considering the coefficients of the body values in the sums 
S [IPj^'PJ, S [IPa;S)Tj]5 it is easy to show that these sums are equal if, and only if, the reciprocal 
relation (4) holds. , And if they are equal, then X^^g [IP^P^] = X/S [IPj^Pj], so that SIPjP^ = when Xj2 is 
not equal to Xj^^, But the reciprocal relation (4) is necessary to the existence of V, 
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for every pair kL This is the condition that the body equations can be derived 
from a single function as in (2). The numbers ii, 4? ••• K ^^^J be regarded as the 
body values of a function of position. It will be denoted by I and called the 
" unifying factor/' because it allows all the body equations to be expressed in terms 
of a single function V. 

For the sake of a certain transformation/ which will be introduced later, V must be 
quadratic and homogeneous, and therefore both body and boundary equations must 
be linear and homogeneous. 

Thus let 

^^^k =^ (^kl^h + Gk2^h + • . . (^kn'^n 



> 



+ ^/d/3i + ej,2^2 + . . . 6/^/3, 



(5). 



And let the boundary equations be 

where the (c)'s, (6)'s, and (/')'s are constants, most of which will commonly be zero. 
In many physical problems we have the integral <^,, given on the boundary. This 
corresponds to a term^o, independent of the (i//)'s, in each of the equations represented 
by (6). However, it is not <^^ which we wish to expand in the form SA/,P^., but the 
difference ((^.h—^m) between the integral and an arbitrary function satisfying the same 
boundary conditions. In ^,,i— <^„ "th^ Xo terms cancel. Thus we are only here concerned 
with functions having boundary equations which are linear and homogeneous in the (i//)'s. 
Next, equation (4) leads to the following relations between the coefficients of (5) 
and (6) : — 

j — s j "= s 

H^ik'^H 2 <^ijjjk = '^k^ki'^hi 2 (^kjjji ' ' ' ' ' • • (' ) 

for every pair k and I, making at most ^n (n — l) equations. Now it is customary in 
physical mathematics to treat body and boundary conditions separately. Let us 
adopt the same course here and treat only the case in which (7) splits into two parts, 

namely, 

"^i^ik ^ '^/A-; ........... (oj 

to be satisfied within the body region, and 

j — s J — •■? 

^Z ^ ^Ijjjk ^^ % ^ ^kjjjl ' • • » » • * ° • \^) 

j=i i=i 

to be satisfied on the boundary. 

In the centre of a sufficiently large table (7) reduces to (8) even in the most general 
case, for there the (e)'s and (/)'s vanish. 

Stated in words, iiC^, = i^fi^ means that the body equations must be able to be 
brought into forms such that if each pair of points in the table be taken in turn, then 
the integral ^ at the first point of the pair, enters into the body equation 5)'<^ at the 



BY FINITE DIFFERENCES OF PHYSICAL PEOBLEMS, ETC. 353 

second, with the same coefficient that the value of <j& at the second enters into the 
body equation at the first. For example, if xjji, t//25 • • • ^n are spaced in order at 
equal intervals along the X axis, and if ^ = d^/dx^^ then 2)^i/;;, = \jjk-i — 2\jjk + \jjk+i 
and (7) is satisfied when the (^ys are simply unity. But if 2) = d/dx, then 
^^\jj^ = /xSi//;;; = xlJk-i—'^k+u SO that Cki is ±1, while C/;, = 0, so it is impossible to find 
a unifying factor. 

Next, as to {{tei^f^j, = ij^e^fji. A particular condition in which these sums become 
identical is if f^i = e^iiaj where a^ is an arbitrary number, independent of the (i//)'s. 

By (6) this equation implies that 

^3 = aj [^i/ii/zi + 62/2^2 + . . . + enjin^n~] (10)- 

Equation (10) states that the boundary values will satisfy their part of (4) if every 
boundary number is formed as follows : Take each body -point. Multiply the value 
of (f) there by the unifying factor and by the coefficient of the said boundary number 
in 2)'^ at this point. Sum this product for all body points and multiply the total by a 
number independent of the body values. This is not the most general way of 
satisfying (9), but it is a common one in physical problems. Thus, for example, if 

the equation to be solved is (^-1 + ^--2)<^« = 0, a boundary condition of frequent 
occurrence is to have ^„ given, and therefore <^r)i—^u = 0. This corresponds to aj = 0. 

Another common condition is = -—^ = -Ax^^ — ^ ^ ^here dn is an element of the 

an on 

normal to the boundary. Then a boundary number will either be equal to a body- 
number or to a weighted mean of two neighbouring body-numbers ; in either case 
af^ — 61/1 + 62/2+ '*'+enjin, in which all but one or two of the (e)'s will vanish. 

Lastly, it will be necessary in what follows that V should be of one sign for all 
values of 1//1, 1//2) • . . ^n- This will be so if 3^V/3i/// is of the same sign for every i//. 
Whence, by (3), (5), and (6), ik{6kk + ^6^jfj^} must be of the same sign for 
every >^ (H)- 

As before, we will treat only the case in which the boundary and body equations 
satisfy this separately, so that every i^fi^^, and every ij^6j,jfjj, have all the same sign. 

For simplicity we will take V as positive (12). 

Then % has the same sign as Cj^j,. Further, as the sign of the operator 2)' does not 
afiect the solution of the equations 2)'i/^i = O...T)'i/^^ = we may alter it arbitrarily, and 
we will for simplicity suppose it chosen at each point so as to make c^k positive . (13). 

The 4... {„ will then also be positive (14). 

For convenience, a table is added giving the forms of I and Y for some common 
forms of 2)'. These were taken by analogy from infinitesimals, and afterwards 
verified by trial. Thus for 2) = V^ we have Y = |- 1|| [V^]^ dr, where V is the vector 
operator Nabla and dr an element of volume. On expressing Y in orthogonal curvilinear 
co-ordinates qiq2<][^ and varying ^, it is found that 8Y ^ \\\ Scj) . V^<^ .1 .dr + a surface 
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integral^ where the unifying factor I is equal to minus the ratio of the volume of the 
elementary co-ordinate block to the product dqidq2dq^,^ 

(3 denotes a summation for every difference which involves a body-number directly 
(that is not by way of boundary equations (6)). 

In the following table the co-ordinate difference is supposed constant for each 
co-ordinate separately : — - 



2)' with sign to make 

^kk plus 






,«x^ ^f "Sz- 









-f- A 



+ 



"& 



^x^ ^x'^y' % 



Unifying 
factor I. 



Unity 



r 



Unity 



Unity 



V. 



(m\<sm\@(y^ 



K'Sx 



%. 



«z 



'■SA 



21 r /^^X2- 



Siri cr- > + SW\- 



'Sz) 



+ 






*@ 



J rp 






At. . 



4 4 4 



^x' %' «^^ 



1 

2 






2 



^+-i^Y 



,^x^ %7 



Note that the averager jx occurs in 2)' but not in V. Having now set out the 
conditions of existence of V in a form in which they can easily be applied to test any 
operator 2)' given with unifying factor and boundary conditions, let us pass on to 
deduce the properties of the principal modes of vibration from the existence ofY, 

Let 



then 



2T = %t|fi^ + %t/r/+...4t|r^^ 



dT/dxjJi = ii^i . 
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(16), 



and by (14) T is essentially positive. 

Now V and T being real quadratic forms in the same variables, one of which, T, is 



* Simply putting "S for 3 in the infinitesimal V does not give the finite-difference V except in special 

* 

cases. — Apil^ 1910, 
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definite,'^ may be transformedf to real sums of squares of new variables A which are 
real linear functions of i/fi..,t|;^, 

2 V = Xi Aj + A2 A2 + . . . + Xyj Ayj. , , (18), 

X/...X/ being the values of X which cause the vanishing of the determinant whose 
(^, il)th element is 

g^P^(V-X^T) , . . . (19). 

But by (3) and (16) this element is 

4.^(3)%.-- XV.). 

And the vanishing of the determinant is, therefore, the condition that the n body 
equations 

2)'(i//i, i//2, ... t/r,) = \'{xjj,, xIj2, ... xjj,) (20) 

should have an integral other than zero. The integrals thus defined as to their body 
values by (20) have already been denoted by Pj, Pg, ... P^. 

The determinant being of the nth degree in X^ vanishes for n values of X^. Since 
one of the forms, T, is definite these n roots are all real (Krofecker, loc. cit), and 
since both T and V are positive Xl^..X/ are all positive ........ (21). 

Now let the coefficients in the transformation of T and V be the (p)'s defined by 



Ai = h\jJiPu + i2\lj2Pl2+ ^'-iiifnP 
As = h\ljlP21 + h^2P22+ '"in^nP 



In 



> 



nn J 



K = k^lPnX + H'^2Pn2 + • • ^^J^nP 

Then differentiating (17) and (18) by i//^ and using (22) 



• • • • 



(22). 



1 7fV 

'l>l = J^^^ = -^lPu + ^2P^l+...+KPnl (23), 

^''/'2= -3— = VAiPu + X/A2^2,+ ...X/A„|)„; . . . . . (24), 

for ^ = 1, 2 ... w. Now if we limit the hitherto arbitrary (t/»)'s by making all the 
(A)'s vanish except Ay, then 

y^ii = AjPji, ^'xjfi = X/AjPji ..... (25), (26), 



"^ ii 



Definite " here means : one-signed and vanishing only when 1^], ^2, ... i'n all vanish. 
t Kronecker, quoted by Bromwich, 'Cambridge Tracts in Mathematics/ No. 3, §26; also by 
Webster, * Dynamics/ Appendix V. 

2 2: 2 
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so that under this limitation + 2)'i//z = A./^^ Now if hjj, is defined as related to pji in 
the same way as yS^. is related to i//^, that is by equation (6) then operating on both 
sides of (25) with 2)^ we have 5)'i//; = Aj^^pji. Substituting in (26) 

^%i = h'PM ■ ■ • (2^)- 

Thus it appears that the (p)'s and (6)'s, as defined by (22) and by (6), are the 
body and boundary values of the integrals Pi ... P^. 

Equation (23) is therefore equivalent to the statement that an arbitrary function 
/may be expanded in the series 

/= AiPi + A,P,+ ...A,P, ........ (28) 

as to its body points, and as / and the (P)'s must all satisfy the same boundary 
condition the same expansion holds good on the boundary also. Analogous to this in 
infinitesimals, and for a special form of 2)', are the expansions of arbitrary functions 
in series of sines, Bessel functions, spherical harmonies, and other integrals of 
(V^+F)(^=0. 

Equations (22) mean that to determine the coefficient of any principal mode of 
vibration in the expansion of an arbitrary function we must multiply the arbitrary 
function by this mode of vibration and by the unifying factor, and add up the products 
at the body points only. This is analogous to the well-known Fourier method of 
determining coefficients. 

Squaring equation (23), multiplying by ii and summing over the body points and 
using (15) and (17), we find 



l=n 1=71 -^ 

t iip\i = 1 , t iiPjiPki = . . . . 

1=1 1=1 



U29) and (30). 



which mean the same as SIP/ = 1, SIP/,P^ = . 



The above proof holds good even if any of X^ .,. X/ are equal to one another. 
This completes the properties of the (P)'s which we require.'^ It remains to consider 
some approximations to Xi^...X/. 
Since by (17) and (18) 

V _ X/A/ + X2^A2^+ ...X^ A^ /q-j\ 



X A^x "4" XX2 "!"••• -^ 



2 
n 



* The transformation of {f)'% to (A)'s resembles an orthogonal transformation of co-ordinates in that 
there are two other relations similar to (29) and (30). These are found by squaring (22) and using (17) 



j = % j = -?! 



and (15). They are 4 2 'p^i? = 1 and 2 ^ji^j2 = 0. The first means that the sum of the squares of the 

values of the n harmonics at the ^th body point is 1/4. The second means that if we take each principal 
mode of vibration, form the product of its values at two fixed body points, and then add up these products 
for all the modes, their sum is zero. — Apil^ 1910. 



BY FINITE DIFFEEENCES OF PHYSICAL PEOBLEMS, ETC. 357 

it follows that when f=zVj, then V/T = X/, and when / differs slightly from P;;. then 
V/T will differ from X/ by a small quantity of the second order (Rayleigh, ' Theory 
of Sound/ § 88). In § 3 '2 '3 a guess is made at P;;., and X/ is then calculated roughly 
as V/T. 

The special values Xi^...X/ for any region are included between the greatest and 
least of the X^ pertaining to any region, which includes the former region when the 
boundary values of both regions vanish. For let P/„ the integral of (2)' + X/) P;;. = 
for the smaller region, be set round about with noughts until a larger boundary is 
reached. The conditions for the existence of a function V will then be satisfied for 
the larger region. Let quantities belonging to the larger region be distinguished by 
dashes. Then by (28) P^ may be expanded in the form 

P;;. = SB^P';, whcre the (B)'s are constants. 

Now by (31) X/ = V/T, when the body values have the values given by (^ = A;;.P;;.. 

Also, by (15), T' = T since the added squares are zero. Also, by (3), ikS^'y\fk is equal 
to both dY/dxfjf, and dY^/dxjj^, so that V and V can only differ by terms independent of 
^ij ^2) • • • ^n- But there are no other variables in Y\ Therefore V = V and X/ = V'/T^ 

Now, expanding V and T' by (17) and (18), 



^k 



\ f2r\ 2 I \ f2r\ 2 I I \ f2r\ 2 

Ax VTx + A2 VT2 + . . . + Ay VjTy 



And therefore X/ lies between the greatest and least of the X'^ If, however, the 
boundary values do not vanish then dY^/Bxjj^ = +4*^'i///, when the (^)'s are regarded as 
independent of the (i//)'s during the differentiation, whereas dY/difj^ = the same 
quantity when V is expressed as a function of the (xfj)'^ entirely, and a more detailed 
examination shows that X/ need not lie between the greatest and least values of X'^ 
for an enclosing boundary with zero boundary values. However, X/ is necessarily 
finite, for it is the root of a rational integral function with fi.nite coefficients. 



